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ABSTRACT 


> 

Presented  is  a  new  method  of  separating  the  zeros  of  a  Finite 
Impulse  Response  (FIR)  filter  producing  an  optimal  digital  filter 
or  surface  acoustic  wave  (SAW)  design  implementation.  Overviews  of 
zero  extraction  algorithms  and  of  FIR  filter  design  using  the  Remez 
Exchange  algorithm  are  presented  (McClellan  et  al .  1973). 

The  computer  aided  design  (CAD)  procedure  presented  allows  the 
designer  to  specify  the  general  filter  characteristic  which  the 
Remez  algorithm  translates  to  FIR  time  domain  coefficients.  These 
coefficients  are  readily  translated  to  the  frequency  (z)  domain, 
producing  an  Nth  order  polynomial  in  z.  The  characteristic 
Dolynomial  is  factored  to  determine  all  roots  or  zeros  using  a 
three-stage  factoring  program  presented  by  M.A.  Jenkins  (1975).  The 
roots  are  optimally  separated  into  two  groups,  each  of  which  is 
recombined  to  form  mutually  exclusive  functions.  The  two 


functions  are  then  implemented  as  transducers  of  a  SAW  device  or 
as  a  two-processor  digital  filter.  The  concept  may  be  extended 


to  more  than  two  subgroups  for  mul ti -processor  digital  filter  designs. 
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CHAPTER  I 
INTRODUCTION 

Numerous  techniques  exist  for  designing  and  implementing 
finite  impulse  response  (FIR)  filters.  Many  of  these  techniques 
can  be  traced  to  antenna  array  design  methods  popularized  in  the 
1940s  and  1950s  (Balanis  1982).  Among  the  more  prominent  antenna 
array  design  techniques  are  the  methods  by  Fourier  transform, 
Schelkunoff  polynomial,  Dolph-Chebyshev,  Taylor  line-source 
(Chebyshev  Error)  and  Woodward.  These  have  spawned  many  of  the 
popular  contemporary  FIR  design  techniques  such  as  the  Remez 
exchange  algorithm  based  on  the  Chebyshev  Error  method,  popularized 
by  McClellan,  Parks  and  Rabiner  (1973)  and  non-iterative 
eigenfunction  synthesis  design,  introduced  by  Devries  (1973)  and 
similar  to  Woodward's  method.  Another  FIR  design  approach  employs 
a  technique  known  as  linear  programming  (Rabiner  1972a ,b). 

Each  of  these  techniques  will  yield  FIR  transfer  functions 
which  may  be  readily  implemented  using  a  surface  acoustic  wave  (SAW) 
device  or  a  digital  filter.  The  SAW  filter,  a  two-transducer 
device,  requires  that  the  transfer  function  be  split  in  some  fashion 
between  the  transducers.  The  digital  filter  may  be  optionally 
implemented  using  two  (or  more)  processors  to  increase  throughput 


rate. 
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The  conventional  approach  to  implementing  the  SAW  device  is, 
basically,  to  construct  one  transducer  such  that  it  contains  the 
entire  FIR  and  to  construct  the  other  transducer  such  that  it 
emulates  a  rect  function.  This  imposes  a  requirement  upon  the  first 
transducer  that  it  be  capable  of  handling  all  of  the  dynamics 
associated  with  the  transfer  function.  Similarly,  a  single-processor 
digital  filter  implementation  demands  that  the  processor  be  able  to 
handle  a  wide  range  of  FIR  coefficients,  as  well  as  all  of  them  at 
once.  These  are  not  necessarily  optimal  implementations. 

Some  efforts  to  evenly  split  the  transfer  function  between  the 
two  transducers  of  a  SAW  filter  have  been  made  by  Morimoto,  Kobayashi 
and  Hibino  (1980)  and  in  work  by  Ruppel ,  Ehrmann-Falkenau,  Stocker 
and  Mader  (1984,  1985).  In  both  cases,  these  teams  split  the 
transfer  function  into  groups  of  alternating  zeros  or  roots  of  the 
transfer  function  about  the  unit  circle.  Any  attempts  to  further 
optimize  the  filters  were  made  at  the  expense  of  altering  the  overall 
frequency  response  in  a  process  called  compensation. 

This  thesis  presents  an  approach  to  near  optimally  split  the 
transfer  function  between  the  two  transducers  or  processors  without 
altering  the  overall  frequency  response.  The  approach  seeks  to 
minimize  non-linear  and  finite  wordlength  error  effects  by  reducing 
the  required  tap  range  for  each  transducer,  or  the  required  range 
of  coefficients  used  in  a  fixed-point  digital  processor. 


CHAPTER  II 

OBJECTIVE  OF  PROPOSED  WORK 


Optima]  FIR  Implementation 
Via  Two-Transducer  Design 


Generally,  the  word  "optimal”  implies  having  attained  a  most 
favorable  condition  or  degree.  Many  parameters  must  be  considered 
during  the  design  of  a  SAW  or  digital  filter.  Addressed  in  this  thesis 
are  those  concerned  primarily  with  filter  order  and  coefficient  dynamic 
range. 

FIR  Coefficients  Via 
the  Remez  Algorithm 

The  first  stage  of  the  design  is  accomplished  using  the  Remez 
exchange  algorithm  (Remez  1957)  to  generate  a  "best  fit"  Chebyshev 
polynomial  to  a  set  of  frequency  response  specifications.  The 
approximation,  and  subsequent  conversion  to  an  impulse  response,  is 
accomplished  by  the  modified  McClellan,  Parks  and  Rabiner  (1973) 
program  presented  in  Appendix  A.  The  program  has  been  altered  to 
permit  the  design  of  filters  of  up  to  an  order  of  one-thousand.  An 
initial  guess  of  optimal  filter  order  is  obtained  using  a  formulation 
presented  by  Vaidyanathan  (1985),  which  states: 


-10  iog^Q  62  -  13 


14.6Af 


(2.1) 


a)  =  stop-band  edge  frequency 


uip  =  pass-band  edge  frequency 


=  pass-band  tolerances 


62  =  stop-band  tolerances 


This  Ng  provides  a  starting  point  for  the  filter  order  in  the  program 


The  program  then  iterates,  increasing  N  each  time,  until  the 
specifications  are  met. 

Once  the  FIR  coefficients  are  found  by  the  modified  McClellan, 
Parks  and  Rabiner  (1973)  program,  they  may  be  easily  arranged  as  the 
coefficients  of  a  z-domain  polynomial,  i.e.,  the  z-transform  of: 


h(n)  =  h(t  -  nT)  =  { 


n  =  0,  1,  2,  ...»  N 


otherwise 


(2.2) 


is 


N 


H(z)  *  Z  h(n)z 
n*0 


-n 


(2.3) 


where: 

ap  =  the  coefficients  generated  by  the  program 

At  this  point,  the  impulse  response  could  be  implemented  as  a  single 
SAW  transducer  or  as  a  single  processor  digital  filter.  In  order  to 
split  up  the  response  between  two  transducers  (or  processors),  H(z) 
can  be  expressed  as: 


5 


v'N-kb 

H(z)  =  - Hjfz)  H2(z) 


(2.4) 


This  equation  shows  all  of  the  poles  of  a  finite  impulse  response  filter 
to  be  at  2  =  0.  All  of  the  filter's  zeros  may  be  found  by  factoring 
the  numerator.  A  judicious  separation  of  the  zeros  (and  poles)  can 
then  be  assigned  to  Hj(z)  and  HgU),  the  two  transducers  of  the  SAW 
filter.  In  a  similar  fashion,  the  transfer  function  could  be  split 
into  Hj  N(z)  for  a  DSP  filter  to  increase  speed  and  dynamic  range. 

Obtaining  the  Zeros  of  the  Characteristic 

This  is  the  most  difficult  phase  of  the  optimization.  The 
factoring  of  high  order  polynomials  was  the  subject  of  considerable 
effort  by  mathematicians  during  the  mid-seventies.  Of  the  many 
techniques  surveyed,  the  Jenkins  and  Traub  (1970)  algorithm,  and 
program  by  Jenkins  (1975),  appeared  to  be  the  best  choice.  This 
algorithm  employs  a  three-stage  process  to  determine  the  roots  of 
an  Nth  order  real  polynomial.  It  is  globally  convergent  and  does 
so  very  rapidly.  The  program  is  extremely  well  written  and  is 
quite  elaborate,  to  the  point  of  compensating  for  specific  machine 
accuracy  limitations. 

The  Jenkins  (1975)  program  factors  the  H(z)  numerator  polynomial 
and  returns  the  real  and  imaginary  portions  of  the  roots  of  H(z). 

The  complex  roots  will  always  appear  in  one  of  two  possible  ways. 

A  set  of  roots  may  appear  as  complex  conjugate  pairs  (quadratic 
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factors),  each  with  a  magnitude  of  one  corresponding  to  the  |z|  =  1 
unit  circle.  These  roots  will  always  correspond  to  the  stopband 
zeros  of  a  filter  design.  Another  form  in  which  they  may  appear  is 
as  a  set  of  two  complex  conjugate  pairs  arranged  symmetrically  about 
the  unit  circle  so  as  to  satisfy  the  condition  that: 

zx  z2*  =  1  (2.5) 


These  zeros  correspond  to  the  passband  zeros  of  the  filter.  A  diagram 
best  illustrates  these  concepts  (see  Figure  1).  Real  roots  may  occur 
on  the  unit  circle  or  in  a  manner  similar  to  the  passband  zero  case. 
Summarizing  the  above,  zeros  of  H(z)  will  occur  as  first,  second  and 
fourth  order  factors. 


Selection  of  Zeros  and  Subsequent  Reconstitution 
Once  the  zeros  of  the  transfer  function  have  been  determined, 
they  must  be  separated  into  sub-groups  and  recombined.  A  simple 
algorithm  which  generates  all  possible  combinations  of  N  roots  taken 
K  at  a  time  is  used  to  separate  the  roots  and  form  the  sub-groups. 

A  polynomial  is  constructed  from  each  group  by  multiplying  the  roots 
within  its  group  and  the  split  design  is  rated  as  to  the 
desirability  of  the  design. 


The  criteria  used  in  the  selection  process  employed  here  seeks 
to  maximize  the  average  tap  or  coefficient  value,  while  minimizing 
the  range  and  variance  of  those  same  values.  These  qualities  are 


v  .  *.  . 
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desirable  in  SAW  devices  since  we  usually  desire  a  maximum  finger 
overlap  and  want  to  avoid,  as  much  as  possible,  large  numbers  of 
very  small  tap  weights  which  may  increase  diffraction  effects.  In 
the  case  of  digital  filters,  these  problems  translate  to  finite  word 
length  problems  and  device  dynamic  range. 

In  order  to  evaluate  the  relative  merits  of  one  combination 
over  another,  the  following  figure  of  merit  is  proposed: 


Design  FOM  = 


R*°x 


(2-6) 


where: 


=  average  of  the  tap  weights  of  both  transducers  combined 
=  the  range  of  the  coefficients 
=  the  variance  of  the  coefficients 


The  ratio  yields  a  figure  of  merit  used  to  rate  a  given  design. 
Obviously,  this  concept  could  be  extended  to  more  than  two 
sub-transfer  functions  for  the  digital  filter  case. 

This  thesis  proposes  to  study  low  order  filters  using  this 
separation/reconstitution  technique  and  to  apply  detected  trends,  if 
any,  in  a  general  sense. 


IVT 
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CHAPTER  III 
FIR  FILTER  DESIGN 

An  excellent  review  of  finite  impulse  response  filter  theory 
is  presented  by  Lawrence  R.  Rabiner  and  Bernard  Gold  (1975)  and  by 
Rabiner,  McClellan  and  Parks  (1975).  This  review  provides  a 
theoretical  background  for  implementing  the  Weighted  Chebyshev 
Approximation  filter  design  technique  via  the  Remez  Exchange 
Algorithm  (Remez  1957).  That  presentation  draws  upon  the  work  of 
Parks  and  McClellan  (1973),  who  devised  a  general  computer  program 
incorporating  the  above.  Their  program  was  used  in  this  thesis 
to  provide  the  FIR  transfer  functions.  An  overview  of  the  theory 
leading  from  FIR  theory  to  a  brief  program  description  is 
presented.  The  following  is  adopted  from  Rabiner  and  Gold  (1975). 

FIR  Filter  Frequency  Response  Overview 
A  finite  impulse  response  describes  a  system  which  can  be 
modeled  by  a  difference  equation  in  the  form: 

M  b 

y(n)  =  E  (r— )  x(n  -  r)  (3. 

r=0  ao 


Since  the  system  output  is  the  convolution  of  the  system  input, 
x(n),  with  the  system  impulse  response,  h(n),  the  impulse  response 
can  be  readily  seen  to  be: 
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n  =  0,  1,  2,  ...»  M 

(3.2) 

otherwise 

The  above  equation  describes  the  discrete  time  domain  coefficients  of 
the  system  FIR.  This  same  response  may  be  described  in  the  frequency 
domain  by  taking  the  Fourier  transform  of  h(n)  to  obtain  H(eJ<i)): 

oo 

H(ejw)  =  Z  h(k)  e~Ja,k  (3.3) 

k=-00 


Since  h(n)  is  finite  in  length  with  respect  to  time,  H(eJa))  must  be 
infinite  with  respect  to  frequency.  However,  for  discrete,  sampled 
systems,  H(eJU))  is  periodic  with  respect  to  the  sampling  frequency, 
i .e. , : 

H(eJU))  =  H[eJ^+27rm^]  m  =  0,  +1,  +2,  ...  (3.4) 


which  is  periodic  in  frequency  with  a  period  of  2ir.  This  fact 
allows  us  to  restrict  our  requirement  to  define  H(eJU))  in  practical 
filtering  applications  in  terms  of  the  sampling  frequency,  consisting 
of  N  samples  over  a  period  equaling  the  length  of  the  time  domain 
Impulse  response  (without  any  augmenting  zeros).  It  also  allows  us 


to  state  that: 


H(eJW)  =  l  h(k)  e'juk 
k*0 


.VAVwV 


■»*  V  *  . 

k*  *  *  ' 


(3.5) 


The  function,  H(eJ  ),  can  be  described  in  terms  of  its  magnitude 
and  phase  as: 


H(eja))  =  +  |H(eja,)|  ej0(a,)  (3.6) 

or 

H(eJa))  =  H  (ej“)  eje(“^  (3.7) 

where : 

H  (e^)  *  a  real  function 

0(w)  *  constrained  to  describe  a  linear  phase  characteristic, 

i.e. ,: 

0(w)  =  -  aw  -ir<_a)<TT  (3.8) 

with  constant  phase  delay  implied  by  the  constant,  -  a.  The  function 
can  be  written  in  trigonometric  form  as: 

+  |H(e0u)|  {cos  (aw)  -  j  sin  (otto)}  (3.9) 

In  order  to  find  a,  equate  the  real  and  imaginary  parts.  Then,  we 
may  describe  cos  (aw)  and  sin  (aw)  as: 


+  ]H(e^u)|  cos  (aw)  = 


N-l 

I  h(n)  cos  (wn) 
n=0 


(3.10) 


and 
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+  |H(eJW)|  sin  (aw)  =  z  h(n)  sin  (cun ) 


N-l 


n=0 


(3.11) 


and  set  up  the  ratio: 


sin  (gw) 
cos_ (c ttii) 


N-l 

Z 

tan  (aw)  *  jjry 
Z 

n=0 


h(n)  sin  (on) 


h(n)  cos  (on) 


(3.12) 


Cross  multiplying: 


N-l  N-l 

Z  h(n)  sin  (aw)  cos  (nw)  -  Z  h(n)  cos  (aw)  sin  (nw)  =  0 
n=0  n=0 

(3.13) 


Using  the  trig  identity: 


sin  (u-v)  =  (sin  u)(cos  v)  -  (cos  u)(sin  v)  (3.14) 


yields: 

N-l 

I  h(n)  sin[(a-n)u]=  0  (3.15) 

n=0 

Equation  (3.15)  is  in  the  form  of  a  Fourier  series.  For  equation 
(3.15)  to  be  valid  for  an  odd  symmetrical  series  (see  Figure  2), 
a  and  h{n)  must  be: 


(3.16) 


h(n)  =  h(N-l-n)  0  <  n  <  N-l  (3.17) 

In  the  case  of  an  even,  symmetrical  series  (see  Figure  3),  a  will  not 
be  an  integer.  The  use  of  the  fractional  delay  obtained  here  is  of 
primary  significance  when  designing  differentiators  and  Hilbert 
transformers.  These  are  not  discussed  here,  but  the  reader  is 
referred  to  Rabiner  and  Gold  (1975)  for  an  in-depth  discussion. 

These  values  for  a  and  h(n)  hold  for  constant  group  delay  and  constant 
phase  filters.  If  constant  phase  delay  (phase  divided  by  frequency) 
is  not  required,  i .e. , : 


H(ejw)  =  ±  |H(ejw)[  e 


-  a)  = 


-  t  |H(eJU))  |  e 


jw»|  J(0-C*o) 


(3.18) 


where  the  phase  delay  is  given  by  -  ^  +  a,  then  similar  development 
(Rabiner  and  Gold  1975)  for  an  odd,  anti -symmetri c  series  (see 
Figure  4)  will  lead  to  the  result  that: 


} 

t 


a 


-  N-l 

-  y 


(3.19a) 


i 


and 


(3.19b) 


h(n)  =  -h(N-l-n) 


0  <  n  <  N-l 


(3.19c) 


Series. 


Again,  the  case  of  an  even,  anti -symmetric  series  (see  Figure  5)  is 
of  primary  interest  in  designing  differentiators  and  Hilbert 
transformers.  Equations  (3.16)  through  (3.19)  suggest  four  general 
classes  that  might  characterize  a  linear  phase  finite  impulse 
response  filter: 

1.  Symmetrical  impulse  response,  N  odd 

2.  Symmetrical  impulse  response,  N  even 

3.  Anti -symmetrical  impulse  response,  N  odd 

4.  Anti -symmetrical  impulse  response,  N  even 

It  is  now  possible  to  describe  H(ejw)  to  account  for  these 
possibilities  in  the  general  relationship: 

H(eja))  =  H(ejw)  ej(B_aw)  (3.20) 

Rabiner  and  Gold  (1975)  next  develop  equations  to  define  H(ejLJ) 
in  terms  of  H(eJW)  for  each  of  the  above  cases.  The  results  of  these 
developments  are  summarized  as  follows: 

Case  1:  Symmetrical  impulse  response,  N  odd 
-  iu)  (N-D/2 

H(eJU)  =  I  a(n)  cos  (con)  (3.21) 

n=0 


with 


a(0)  =  h[(N-l)/2] ,  and 
a(n)  =  2h[^-^-  -  n]  for  n  =  1,  2,  ..., 


(N-l)/2 
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which  yields: 


H(eJ") 


e-j«(N-l)/2  »-*>■ '2  a(n)  cos  U) 

n=0 


Case  2:  Symmetrical  impulse  response,  N  even 


-  in  N/2 

H(eJW)  =  Z  b(n)  cos  Mn-0.5)] 
n=l 


with 


b(n)  =  2h(N/2  -  n),  n  =  1,  2 . N/2 


H(eJU))  =  e“ju{N“^/2  V  b(n)  cos  Mn-0.5)] 

n=l 


Case  3:  Anti -symmetrical  impulse  response,  N  odd 


-  io,  (N-D/2 
H(eJU)  =  Z  c(n)  sin  (con) 
n=l 


wi  th 


c(n)  =  2h[(N-l)/2  -  n],  n  =  1,  2,  ....  (N-l)/2 


H(e^u)  «  e-Jw(N-l)/2  e^/2  ^  z*  c(n)  sin  (con) 


n=l 


(3.23) 

(3.24) 

(3.25) 

(3.26) 

(3.27) 

(3.28) 

(3.29) 


Case  4:  Anti -symmetrical  impulse  response,  N  even 


~  s  N/2 

H(eJaj)  =  z  d(n)  sin  [o>(n-0.5)]  (3.30) 

n=l 


with 

d(n)  =  2h(N/2  -  n),  n  =  1,  2 . N/2  (3.31) 

N/2 

H(eju)  =  ej7r/2  z  d(n)  sin  [w(n-0.5)]  (3.32) 

n=l 

Weighted  Chebyshev  Approximation 
The  frequency  response  of  the  desired  system  is  completely 
described  by  H(eJOJ) .  From  the  development  in  the  first  section  of 
this  chapter,  it  is  obvious  that  this  description  for  each  case  can 
be  considered  as  a  series  of  sine  or  cosine  functions.  These  series 
can  be  easily  related  to  Chebyshev  polynomials. 

The  Chebyshev  polynomial  represents  an  expansion  of  cos  (mu)  for 
any  value  of  m.  We  know  that  any  real  function  can  be  represented 
as  a  sum  of  sinusoids.  These  sinusoids  are  of  the  form  cos  (mu), 
with  m  indicating  the  highest  harmonic  required  to  reconstruct  the 
original  function.  Of  course,  some  functions  require  that  m  approach 
°°.  Within  given  limits,  however,  it  is  possible  to  represent  a 
desired  frequency  response  curve  as  a  sum  of  Chebyshev  polynomials 
of  finite  length  m.  The  Chebyshev  polynomial  expansions  take  the 
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m 

cos  mu 

0 

1 

= 

1 

1 

cos  u 

= 

cos  u 

2 

cos  2u 

= 

2  cos2  u  -  1 

3 

cos  3u 

s 

3 

4  cos  u  -  3  cos 

=  cos 

(u),  or  u  = 

cos"1 

(z),  then: 

m 

cos  mu 

Chebyshev 

Designation 

0 

1 

Vz> 

1 

2 

Tjtz) 

2 

2z2  -  1 

T2(z) 

3 

4z3  -  3z 

T3(z) 

(3.33) 


(3.34) 


A  recursive  relationship  emerges: 


Tm(z)  =  cos  [m  cos-1  (z) ]  =  cos  (mu),  -1  <  z  <  1 


(3.35) 


In  essence,  we  are  using  a  sum  of  Chebyshev  polynomials  to 
curve-fit  to  the  desired  frequency  response  from  zero  to  one-half 
the  sampling  frequency.  Given  the  four  cases  discussed  at  the  end 
of  the  first  section  of  this  chapter,  and  using  Chebyshev  polynomials 
to  represent  the  series,  the  problem  defaults  to  determining  the 


scaling  of  the  coefficients  with  respect  to  frequency.  This  process 
is  called  "weighting"  the  approximation  and  is  discussed  in  depth 
by  Rabiner  and  Gold  (1975).  It  is  reiterated  here  briefly. 

For  the  four  cases  described  in  the  first  section  of  this 
chapter,  a  general  expression  can  be  written  to  define  HCe^)  as: 

HCe^)  =  e“’^N"1)/2  eJU/2)L  H(e0w)  (3.36) 

The  exponent  L  will  take  on  a  value  of  either  0  or  1,  depending  upon 
the  case  considered.  Now,  a  table  can  be  constructed  which  shows 

A  • 

values  for  L  and  the  form  of  H(eJW)  for  the  appropriate  case  of 
symmetry  and  N  (see  Table  1).  The  previous  discussion  of  the  form 
of  the  Chebyshev  polynomial  would  suggest  that  the  expressions  for 

A  , 

H(eJU)  may  be  converted  to  summations  involving  cosines  (as  opposed 
to  sines)  using  ordinary  trigonometric  identities.  Once  this  is 
done.  Table  1  can  be  rewritten  in  terms  of  functions  which  are  fixed 
functions  of  oj,  which  will  be  referred  to  as  Q(eJW),  and  as 
functions  of  the  cosine  series,  which  will  be  referred  to  as  P(eJU) 
(see  Table  2).  For  cases  2  through  4,  Q(eJU>)  is  constrained  to  be 
zero  at  either  oj  a  0  or  w  *  tt,  or  both. 

Now,  it  is  possible  to  set  up  a  relationship  between  the 
desired  response  at  given  frequencies  to  within  a  prescribed  accuracy. 
To  do  so,  let  D(eJa))  represent  the  desired  response  of  the  filter 
and  let  W(eJU))  represent  the  weighting  on  the  allowable  error  as  a 
function  of  frequency  regions  or  bands  (i.e.,  the  ratio  of  the 


TABLE  1 


response,  N  even 
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1 


stopband  ripple  to  the  passband  ripple).  With  these  functions,  the 


error  for  a  given  approximation  can  be  calculated  as: 


E(eju))  =  W(eju))  [D(eju))  -  H(ejw)] 


(3.37) 


A  •  A  • 

with  H(eJU))  being  the  trial  design.  H(ejw)  can  be  separated  into  its 


two  symbolic  parts  to  yield: 


E(eja))  =  W(ejaJ)  {D(ejaJ)  -  P(eja))  Q(ejw)] 


(3.38) 


Q(eJU)  may  be  factored  out  of  the  quantity  in  parentheses  since  it  is 
a  fixed  function  of  frequency.  This  yields: 

E(ej(i))  =  W(eju))  Q(eju))  {D(eju)/Q(ejw)  -  P(ejw)]  (3.39) 
Defining  [W(e^“)  Q(e^“)]  as  W(e^),  and  [D(e^aJ)/Q(eJU))  ]  as  Dfe^), 


equation  (3.39)  may  be  rewritten  as: 


E(ejw)  =  W(ejw)  [D(ejuj)  -  P(eju>)] 


(3.40) 


The  problem  now  defaults  to  finding  the  values  for  the  coefficients 
of  the  Chebyshev  polynomials  £P(eJa)) ]  such  that  the  maximum  error 
over  each  specified  frequency  band  is  minimized.  To  accomplish  this, 
the  Alternation  Theorem  is  used.  It  states  (Rabiner  and  Gold  1975): 

Theorem:  If  P(eJW)  is  a  linear  combination  of  r  cosine 
functions,  i .e.,: 


P(eJa))  =  z  ot(n)  cos  (nco) 
n=0 


(3.41) 


•:  •.'’•.  is-:-/  ' . 
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then  a  necessary  and  sufficient  condition  that  P(eJw)  be  the  unique, 
best  weighted  Chebyshev  approximation  to  a  continuous  function  D(eJU)) 
on  A,  a  compact  subset  of  (0 ,tt)  ,  is  that  the  weighted  error  function 
E(ejW)  exhibit  at  least  (r+1)  extremal  frequencies  in  A;  i.e.,  there 
must  exist  (r+1)  points  in  A  such  that  ancl 

such  that  E(e^Wi)  =  -E(eJWi+i),  i  =  1,  2,  ....  r,  and  |E(eJU)i)|  = 
max  [E(e^)]  for  all  u  in  A. 

Rabiner  and  Gold  (1975)  show  that  for  the  four  cases  of  filter 
design  presented,  that  the  number  of  extremal  frequencies  in  H(eJli)) 
obey  the  following  constraints: 


Case  1: 

Ne  <  (N+l)/2 

Case  2: 

Ne  <  N/2 

Case  3: 

Ne  <  (N-l)/2 

Case  4: 

Ne  <  N/2 

(3.42) 


The  extremal  frequencies,  or  extrema,  are  divided  up  between  the  stop 
and  passbands  of  the  filter  and  they  describe  the  peaks  and  troughs 
of  the  Chebyshev  approximation.  A  diagram  best  describes  the 
relationship  of  the  extremal  frequencies  and  the  shape  of  the 
Chebyshev  approximation  waveform  (see  Figure  6). 

There  are  several  ways  in  which  to  obtain  the  extremal 

A  • 

frequencies  of  H(eJU).  The  first  method,  described  briefly,  was 
originally  proposed  by  Herrmann  and  Schuessler.  The  method 
capitalizes  on  the  fact  that  a  local  maxima  (+6)  or  a  local  minima 


(-6)  occurs  in  the  region  of  an  extrema,  and  that  the  derivative 
is  zero  at  that  point.  Two  equations  in  Nfi  with  two  Ng  unknowns 

A  2 

(Ng  impulse  response  coefficients  and  Ng  frequencies  where  H(eJaj) 
obtains  an  extremal  value)  are: 


Jw-  ,  Jaw 

H(e  - 77; —  +  D(e  ’)  i  -  1,  2 . N  (3.43) 

joj.  e 

W(e  M 


(3.44) 


where  E(eJW>)  *  +  6,  and  these  are  solved  iteratively  for  values  of 


Ng.  This  procedure  works  well  for  filters  with  an  order  about  60 
or  less. 


Another  method  used  is  one  devised  by  Hofstetter,  Oppenheim  and 
Siegel  which  is  called  the  Polynomial  Interpolation  Solution  by 
Rabiner  and  Gold  (1975).  The  basic  idea  behind  this  algorithm  is 
that  an  initial  guess  of  the  extremal  frequencies  is  made  and  H(eju)) 
is  evaluated  at  these  points.  The  algorithm  then  searches  for  the 
actual  extrema  found  during  that  trial  and  iterates  again,  this 
time  using  the  newly  found  extrema.  Eventually,  the  process  converges 
to  the  minimum  ripple  attainable  for  a  given  Nfi.  Very  large  order 
filters  can  be  designed  by  this  method.  The  Polynomial  Interpolation 


Solution  technique  is  very  similar  to  the  last  technique  to  be 
described,  the  Remez  Exchange  Algorithm. 
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The  Remez  Exchange  Algorithm 

We  have  seen  that  the  goal  of  obtaining  the  desired  response 
D(eju)  is  met  by  obtaining  the  approximating  function  P(eJa))  which 
best  minimizes  the  weighted  error  function  E(eJa)).  The  Remez 
Exchange  Algorithm  accomplishes  this  by  using  a  dense  grid  of 
frequency  points  to  find  the  extremal  frequencies.  An  initial 
guess  as  to  the  location  of  the  (r+1)  frequencies  is  made,  similar 
to  the  Polynomial  Interpolation  Solution  method.  Then,  the  error 
function  is  forced  to  have  a  value  of  +  5.  The  signs  alternate, 
since  the  extrema  are  expected  to  alternate  above  and  below  the 
indicated  level  by  5  in  the  final  design.  These  constraints 
generate  the  separate  error  function  [E(eJu))]  equation  for  each 
extremal  frequency,  given  from  equation  (3.40)  as: 

W(e  K)  [D(e  K)  -  P(e  K)J  =  (-I)K  6  k  =  0,  1 . r  (3.45) 


which  generates  an  (r+1)  x  (r+1)  matrix  of  equations  to  solve. 

Remez  (1957)  found  an  alternative  closed-form  solution  (appropriately 
modified  for  the  current  variable  set)  to  be: 


aQD(e  )  +  a^e  L)  +  ...  +  arD(e  ) 


aQ/W(e  °)  -  3^(6  x)  +  ...  +  (-l)r  ar/W(e  r) 


(3.46) 


where: 


1 

<xk  -  n' 


(3.47) 
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and 


x.j  =  cos  cj1- 


(3.48) 


At  this  point,  the  optimum  6  for  a  given  set  of  extremal 
frequencies  is  known.  The  next  step  is  to  form  the  approximating 
function  P(eJa>)  along  the  r  extrema  points  by  using  the  barycentric 
form  of  the  Lagrange  interpolation  formula: 


P(eju3) 


k=0  k 

rZl  Bk 
k=0  tx  “  xk} 


(3.49) 


where: 


*k 


r-1 

n 

i=0 


(3.50) 


and 


~  JWk  i. 

C.  =  D(e  K)  -  (-l)k 


~  JWl. 

w(e  k) 


k  =  0,  1,  ...»  r-1 


(3.51) 


A.  =  COS  COj 

x^  =  cos  co^  (3.52) 


x  =  cos  co 


Once  the  approximating  function  P(eJa))  has  been  formed,  it  is 
possible  to  evaluate  E(eJaJ)  along  a  dense  set  of  frequencies  which 
are  equally  spaced  along  the  frequency  axis  from  zero  to  one-half 
the  sample  frequency.  If: 

|E(eju))|<6  (3.53 

then  an  optimal  approximation  to  the  desired  frequency  response  has 

been  found.  If  the  weighted  error  function  exceeds  6,  then  a  new 

set  of  (r+1)  extremal  frequencies  is  chosen  by  selecting  the  peaks 

of  the  error  curve.  This  process  quickly  forces  6  to  converge  to 

its  maximum  value  for  a  given  number  of  extremal  frequencies.  If 

there  are  more  than  (r+I)  extrema  in  E(eJW),  then  the  new  number  of 

extrema  is  retained  and  used  in  the  next  iteration  of  the  process. 

The  final  impulse  response  coefficients  are  obtained  by 

performing  a  2M  point  Inverse  Discrete  Fourier  Transform  on  P(eJU)), 

M 

where  2  >_  N.  Note  that  this  N  is  the  filter  order  plus  1. 


CHAPTER  IV 


ZERO  EXTRACTION  TECHNIQUES 

The  problem  of  extracting  the  zeros  of  a  polynomial  turns  out 
to  be  far  from  simple,  sparking  the  interest  of  mathematicians 
and  scientists  for  centuries.  With  the  advent  of  the  digital 
computer,  the  factoring  of  high  order  polynomials  has  become 
possible,  though  not  entirely  without  grief.  Some  of  the  more 
prominent  approaches  and  associated  problems  are  briefly  discussed 
here. 


Polynomial  Theory 

A  polynomial  in  z  is  an  equation  which  takes  the  form: 

VN  +  v/'1  +  +  V2  +  aiz  +  ao  <4-‘> 

or,  alternatively. 


N 

E  a  zn  (4.2) 

n=0  n 

where  the  coefficients  a^,  a^,  ....  aQ  are  real  numbered  constants. 
This  form  of  equation  is  readily  identified  with  the  equation 
describing  the  finite  impulse  response  filter  z-domain  representation. 
This  polynomial  can  be  expressed  also  as  a  product  of  its  roots  or 


zeros  as: 
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N 

n  (z  -  z  )  (4.3) 

n=l 

Notice  that  for  a  polynomial  with  N  roots,  there  are  N  product-form 
terms  and  N+l  summation -form  terms.  This  can  cause  some  confusion 
at  times.  For  example,  the  McClellan,  Parks  and  Rabiner  (1973) 
program  discussed  in  the  previous  chapter  displays  a  filter  order  of 
N  when,  in  fact,  N  coefficients  are  actually  meant. 

Zero  Characteristics  of  FIR  Filters 

When  described  by  a  polynomial  in  the  z  =  eJW  plane,  FIR  filter 
zeros  plotted  in  the  z-plane  always  have  a  distinct  appearance. 

All  of  the  stopband  zeros  will  occur  exactly  on  the  unit  circle 
and  will  always  occur  as  complex  conjugate  pairs,  unless  they  are 
real.  The  complex  passband  zeros  will  always  occur  in  sets  of 
four,  one  inside  the  unit  circle,  another  outside  such  that  the 
magnitude  of  one  multiplied  by  the  other  will  equal  exactly  one. 

This  pair  also  has  corresponding  complex  conjugates,  hence  the  set 
of  four. 

Due  to  the  nature  of  the  Chebyshev  polynomial  approximation 
of  the  FIR  frequency  response,  there  are  no  repeating  zeros  or 
multiple  roots  to  contend  with.  However,  this  does  pose  problems 
in  other  factoring  situations  and  is  discussed  briefly. 

Factoring  Methods 

Several  unique  approaches  to  zero  extraction  exist.  The  first 
was  by  none  other  than  Sir  Issac  Newton  (1642-1727).  Since  that 


time,  other  algorithms  by  Bairstow,  Lin,  Muller  and  Birge-Vieta 
have  arrived  (Ralston  and  Wilf  1960).  These  algorithms  have  the 
relative  liability  of  not  being  able  to  assure  convergence  for  any 
initial  guess  of  a  root.  Other  methods  which  virtually  assure 
convergence  within  a  class  of  problems  are  methods  of  Lehmer, 

Graeffe  and  Bernoulli.  The  Bisection  method  is  probably  the  most 
crude  method  of  root-finding,  relying  upon  a  purely  iterative  process 
of  testing  discrete  points  in  the  z-plane  until  the  roots  are 
found.  The  latter  four  cases  have  the  disadvantage  of  slow 
convergence.  There  exist  several  matrix-oriented  computer  program 
packages,  such  as  EISPACK  (Smith  et  al .  1976),  which  are  designed 
to  find  the  eigenvalues  of  a  matrix,  another  means  of  finding  the 
zeros.  However,  these  programs  are  extremely  memory-inefficient. 

The  method  adopted  for  this  thesis  project  is  the  Jenkins  and  Traub 
(1970)  method  which  combines  many  of  the  above  techniques,  as  well 
as  ones  by  Traub,  into  a  very  complex  algorithm  and  computer  program 
which  is  convergent  for  a  wide  class  of  problems  and  is  very 
machine-efficient.  Discussed  briefly  are  the  more  prominent  methods. 

Newton's  Method 

The  process  of  finding  a  root  by  Newton's  method  is  perhaps 
the  best  known  and  most  easily  understood  (Blomquist  1968).  It 
is  based  on  beginning  with  an  initial  guess  for  the  root  and 
iteratively  converging  to  the  actual  root  with  the  method  of  steepest 
descent.  The  following  relationship  describes  the  method: 


35 


zn+l  =  zn  -  p(zn>/pl(zn>  <4'4> 

where  zn  is  the  current  guess  of  the  root,  P( zn )  is  the  value  of 
the  polynomial  at  zn,  P*(zn)  is  the  value  of  the  derivative  of 
P(zn)  and  zn+1  is  the  next  guess  Cor,  eventually,  the  root).  This 
process  may  be  carried  out  iteratively  to  any  practical,  desired 


accuracy.  The  root  is  obtained  when  the  difference  between  zn+^ 


and  zn  is  less  than  the  required  accuracy.  Figure  7  graphically 


shows  how  successive  iterations  ultimately  converge  to  the  root. 

Problems  with  this  technique  occur  since  it  has  no  direct  way 
of  determining  if  a  multiple  root  exists  and  the  method  does  not 
always  converge  to  the  root.  An  example  of  how  the  method  may 
fail  is  shown  in  Figure  8.  This  case  demonstrates  that  choosing 
an  initial  guess  too  far  from  the  actual  root  may  prevent  convergence. 
In  this  example,  the  method  will  oscillate  between  z^  and  zQ,  since 
each  point  represents  the  other's  successive  approximation.  The 
method  also  requires  the  use  of  complex  arithmetic  in  evaluating 


P(zn)  and  PA(zn)  in  the  case  of  complex  roots,  which  further  limits 


this  method. 


Bairs tow's  Method 

Prior  to  the  Jenkins  and  Traub  (1970)  approach,  the  Bairs tow 
method  was  regarded  as  one  of  the  best  techniques  for  extracting 
zeros.  The  primary  advantage  which  this  method  has  over  Newton's 
Method  is  that  it  uses  only  real  arithmetic  to  evaluate  the 
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Figure  8.  Newton's  Method  -  Failure  to  Converge. 
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polynomial.  The  basic  idea  is  to  use  Newton's  Method  to  find  unique 
quadratic  factors  to  the  polynomial  using  only  real  arithmetic, 
remove  the  quadratic  via  synthetic  division  and  use  the  well-known 
quadratic  formula  to  extract  the  complex  roots  of  the  quadratic 
factor.  This  technique  automatically  locates  multiple  roots,  since 
it  deflates  the  polynomial  by  an  order  of  two  for  each  quadratic 
factor  found. 

Simons,  Weeks  and  Kotick  (1983)  developed  an  elegant  formulation 
which  best  expresses  Bairstow's  Method.  The  algorithm  begins  with 
a  polynomial  in  the  form  of: 


v 


+  a  2N_1  +  a  7N‘2 

+  aN-lZ  aN-2Z 


+  ajZ1  +  aQz°  (4.5) 


Newton's  Method  is  used  to  approach  the  root  by  using  equation  (4.4). 
However,  Bairstow's  Method  evaluates  PN(z)  and  its  derivati ve(s)  at 
the  pair  of  complex  points  zp  and  zn*  by  using  only  real  arithmetic 
as  follows: 


Let  PN(z)  = 


N 

E 

n=0 


Let  the  quadratic  factor  take  the  form: 

z2  +  ctz  +  6 


(4.6) 


(4.7) 


where : 


a  =  -  2a 


Then,  dividing  Pj^(z)  by  this  quadratic  factor  yields  a  polynomial 

of  order  N-2  with  a  remainder  R,z  +  R  ,  i.e., 

i  o 


Pn(2) 


R,z  +  R_ 


Nv  '  1  *•  rYi 

_  p(z)  + 

z  +  az  +  3  z  +  az  + 


z  +  az  +  3 


(4.8) 


Multiplying  both  sides  by  the  quadratic  factor  yields: 


PN(z)  =  +  az  +  2)  +  RjZ  +  RQ 


(4.9) 


Obviously,  if  RjZ  +  RQ  =  0,  then  the  roots  are: 


-  a  -  a  -  43 


(4.10) 


by  the  quadratic  formula.  Therefore,  the  problem  defaults  to  choosing 

values  for  a  and  3  such  that  the  remainder  is  zero.  R,  and  R  can 

l  o 

be  related  to  P|Sj(z )  by  the  following: 


a/-2  *  (aN-i^aN-i)zN~3^[aN_2-«aN.1+(a2+6)aN]2N~3^ 
z2  *  02  *  6  I  VN  +  aN-l2N_1  +  aN-2zN‘2  +  aN-3z"'3  +  ••• 


(aN-j-aaN)zN  '♦(aN.2-6aN)zN  2+aN.3zN  3+... 

(aN-l~aaN^z  ^+^ctaN-i“aaN)2^  2+^aN_j“a^afj)z^  ^ 
IaN_2-aaN_i+(a2-B)aN]zN"2 

+^aN-3~^aN-l+ct^aN^z  ’*’••• 

(4 

Letting: 


bN-2  =  aN 

bN-3  *  aN-l  "  aNa  =  aN-l  ’  bN-2 

(4 

bN-4  =  aN-2  ’  aNa  ‘  aN-la  "  aNa 
=  aN-2  "  bN-2  ~  bN-3a 
a  recursive  relationship  emerges: 

bN-2-n  =  aN-2-n  '  SbN-2-n  "  abN-3-n  (4 


Iterating  to  n  *  N  yields: 


2~  "  -  b  ,z_x  +  b_2z“^  (4.15) 

r-  +  ctz  +  3  1  * 

then 

RjZ  +  RQ  =  (z2  +  az  +3)  (b^z-1  +  b_2z~2)  (4.16) 

ab  0  3b  .  3b  , 

*  b-lZ  +  b_2  +  ab_j  +  — f-  +  — Ji-  +  — p-  (4.17) 

Equating  like  coefficients: 

R!  *  b-l 
Ro  =  b-2  +  ab-l 

(4.18) 

■  aQ  -  6b0  -  ab^  +  ab_j 

■  ao  -  8bo 


The  polynomial  PN(z)  and  its  derivative(s)  may  now  be  evaluated  at 
any  pair  of  complex  conjugate  points  specified  by  the  chosen  quadratic 


factor.  Newton's  Method  now  follows  easily,  unless  P^z)  evaluates 
to  zero.  If  this  is  the  case,  multiple  roots  exist  for  the  chosen 
pair,  o  +_  ju>.  In  this  case: 


fuhz) 


(4.19) 


is  replaced  by  successive  derivatives 


Ai 


(4.20) 


until  PN  1(zn)  f  0.  Now  the  root,  as  well  as  its  multiplicity,  is 


known . 


The  final  step  of  the  process  is  to  remove  each  quadratic  factor 
from  the  polynomial  (polynomial  deflation)  using  synthetic  division. 
Bairs tow's  Method  is  again  applied  to  the  resulting  polynomial  until 
all  of  the  roots  are  found. 

The  disadvantage  of  Bairstow's  Method  lies  in  the  necessity 
to  make  a  reasonably  accurate  guess  of  the  quadratic  factor.  A  bad 
guess  can  prevent  convergence  in  the  same  manner  as  happens  in 
Newton's  Method.  Both  of  the  above  processes  suffer  from  machine 
rounding  problems  which  occur  with  successive  deflation  of  the 
original,  high  ordered  polynomial. 


Jenkins  and  Traub  Method 


This  is  a  highly  complex  method  which  attempts  to  incorporate 
all  of  the  above  advantages  while  avoiding  the  mentioned  pitfalls. 

The  program  incorporating  the  algorithm  was  the  subject  of  Jenkins' 
doctoral  dissertation  (Jenkins  1975).  The  process  incorporates 
three  stages  of  zero  extraction  using  Bairstow's  Method,  Newton's 
Method  and  three  shifting  techniques  used  to  hasten  convergence. 

The  method  assures  rapid  convergence  for  a  wide  class  of  polynomials. 
Zeros  are  removed  in  roughly  increasing  order  of  modulus;  i.e.,  the 
zeros  closest  to  the  origin  are  generally  removed  first,  the  ones 
furthest  from  the  origin  are  removed  last.  This  is  done  in  order 
to  reduce  the  instability  problems  which  may  accompany  the  deflation 
process.  A  discussion  of  the  variable  shift  algorithm  is  beyond  the 
scope  of  this  thesis  and  the  reader  is  referred  to  the  Jenkins  and 
Traub  (1970)  paper  for  a  formal  theoretical  treatment. 

One  of  the  interesting  aspects  addressed  by  Jenkins  in  writing 
the  program  is  that  it  takes  into  account  the  specific  capabilities 
and  limitations  of  floating  point  manipulations  on  a  given  machine. 
This  feature  allows  the  program  to  be  customized  to  a  specific  machine 
in  order  to  achieve  the  highest  possible  zero  extraction  accuracy 
for  that  machine  (using  the  Jenkins  and  Traub  algorithm). 

The  program  in  Appendix  B  appears  to  be  the  current  state-of- 
the-art  in  non-matrix  methods  of  polynomial  factoring.  For  example, 
the  1986  IMSL  math  libraries  make  use  of  this  algorithm  for  their 
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zero  extraction  approach.  Schelin  (1983)  also  indicated  that  this 
was  the  prime  non-matrix  type  algorithm  as  of  1982. 

In  practice,  the  program  handles  up  to  roughly  an  order  of 
one-hundred  with  little  difficulty.  Convergence  problems  begin  to 
occur  with  increasing  frequency  beyond  this  limit,  based  upon 
actual  tests. 


CHAPTER  IV 


ZERO  SEPARATION  AND  RECONSTITUTION 

The  prime  reason  for  judicious  zero  separation  is  based  on  a 
desire  to  increase  the  dynamic  range  of  the  transducer  or  DSP  device 
without  sacrificing  any  of  its  transfer  characteristics.  The 
dynamic  range  is  largely  affected  by  relatively  small  FIR 
coefficients.  These  small  coefficients  correspond  to  small  area 
overlaps  of  fingers  in  SAW  devices  and  to  small  register  coefficients 
in  DSP  filters.  The  net  effect  in  the  SAW  device  is  for  the  small 
overlap  to  appear  more  like  a  point  source  wave  generator  as  opposed 
to  a  desired  planar  wave  source.  Second  order  effects  also  begin  to 
become  more  predominant  for  this  situation.  Similarly,  DSP  filters 
suffer  from  rounding  effects  when  forced  to  sum  products  of  very 
large  and  very  small  filter  coefficients,  even  when  floating-point 
arithmetic  is  used.  The  following  DSP  example  demonstrates  this 
problem. 

1.  Let  the  internal  number  representation  be  from  -  1.00  to 

1.00. 

2.  Let  the  system  function  be: 

y(n)  =  (1.00)  x(n)  +  (0.02)  x(n-l) 
which  is  an  FIR  filter  with  coefficients  1.00  and  0.02. 


46 


3.  Let  x(0)  =  x(l)  =  0.20  and  x(-l)  =  0.0,  then  it  follows 
that: 

y(0)  =  (1.00)  (0.20)  +  (0.02)  (0.00)  =  0.20 

y(l)  =  (1.00)  (0.20)  +  (0.02)  (0.20)  =  0.204  (actual) 

However,  y ( 1 )  =  0.204  will  be  truncated  to  0.20,  since  the  internal 
precision  only  allows  two  decimal  places. 

In  practice,  it  may  not  be  possible  to  completely  avoid  the 
problems  associated  with  dynamic  range,  but  it  is  desirable  to  attain 
the  best  dynamic  range  for  a  given  design.  The  above  example 
demonstrates  that  three  qualities  of  the  FIR  coefficients  bear  close 
scrutiny  during  the  design  process:  the  average,  the  variance  and 
the  range  of  the  coefficient  values. 

Statistical  Qualities 

The  highest  average  value  for  the  FIR  coefficients  provides  the 
greatest  average  finger  overlap  on  a  SAW  device,  thus  the  highest 
average  energy  injection  into  (or  removal  from)  the  substrate. 
Similarly,  the  highest  average  coefficient  value  in  a  DSP  filter 
produces  data  with  the  highest  average  value  within  the  register 
working  ranges.  Since  sign  changes  may  be  handled  easily  in  either 
type  of  device,  the  average  is  taken  of  the  absolute  value  of  the 
FIR  coefficients.  The  average  is  defined  as: 

N 

a  ^  |h(n)|  (5.1) 

tr*  n=l  _ 
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The  variance  of  the  FIR  coefficients  indicates  how  much  they 
change  on  the  average.  In  other  words,  the  variance  indicates  how 
smooth  the  overall  envelope  is.  Again,  since  sign  changes  may  be 
handled  easily  in  either  type  of  device,  the  variance  is  taken  of 
the  absolute  value  of  the  FIR  coefficients.  The  variance  is  defined 
as: 


2 
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N  ,  N 

N  x  Z  Ih(n):r  -  {  I 
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|h(n)|}‘ 


TT(N-T) 


(5.2) 


The  range  of  the  coefficients  is  an  absolute  indication  of  now 
far  apart  the  minimum  and  maximum  coefficient  values  are.  Since  the 
coefficients  are  normally  scaled  to  a  maximum  of  1.0,  the  range  will 
provide  an  indication  of  the  minimum  coefficient  value.  As  in  the 
above  two  statistical  qualities,  the  range  is  taken  for  the  absolute 
coefficient  values.  It  is  defined  as: 

Range  ■  |h(n)max|  -  |h(n)m1„|  (5.3) 

These  three  qualities  allow  the  designer  some  means  of  rating  one 
given  design  against  another  quantitatively,  leading  to  a  design 
Figure  of  Merit  (FOM). 


The  Figure  of  Merit  (FOM) 


A  design  objective  which  requires  the  best  split  of  the  zeros 
of  the  transfer  function  requires  that  the  designer  be  able  to  rate 


one  design  over  another  until  the  best  is  found.  This  can  be  done 
by  assigning  an  FOM  to  the  design  based  on  the  three  statistical 
qualities  discussed  above.  The  goal  addressed  here  is  to  split  the 
zeros  between  two  transducers  or  processors  such  that  a  gain  in 
available  dynamic  range  is  found.  To  do  this,  the  following  procedure 
relating  the  combined  transducer  coefficient  average,  variance  and 
range  values  to  a  figure  of  merit  is  proposed: 

1.  Normalize  all  hj(n)  coefficients  to  1.0  maximum. 

2.  Normalize  all  h2(n)  coefficients  to  1.0  maximum. 

(NOTE:  For  the  special  case  where  all  of  the  FIR 
coefficients  are  implemented  by  h,(n),  then 
h2(n)  is  set  to  1  to  represent  the  impulse 
function  since  the  convolution  of  the  impulse 
response  with  an  impulse  is  the  impulse 
response.) 

3.  Form  an  array  consisting  of  Ih^n))  . 

4.  Concatenate  the  |h2(n)|  values  to  this  array. 

5.  Find  the  average  (x)»  variance  (a2)  and  range  of 
the  newly- formed  array. 

6.  Apply  the  relation: 


This  equation  reflects  a  desire  to  maximize  the  average  while 
minimizing  the  variance  and  range  of  the  coefficient  values.  It 
provides  the  designer  with  a  means  of  rating  one  set  of  split  zeros 


versus  another. 


In  order  to  determine  if  an  optimum  zero-splitting  pattern 
might  exist,  all  possible  combinations  for  distributing  the  zeros 
must  be  made  and  each  combination  tested  using  the  FOM  procedure. 
This  process  is  very  tedious  but  may  be  readily  implemented  on  a 
computer  due  to  the  iterative  nature  of  the  process.  Since 
there  are: 


total,  non-repeating  combinations,  the  process  has  a  practical 
computational  upper  limit  of  about  FIR  order  N  =  40  (which  has  over 
10**  combinations)  on  a  VAX  11-759.  However,  equations  of  low  order 
may  be  used  to  model  any  trends  applicable  to  the  higher  order 
cases  encountered  in  SAW  devices. 


Computer  Implementation 

Appendix  C  shows  a  listing  of  the  program  COMBO  used  to  generate 
(Beckenbach  1964)  and  rate  all  of  the  combinations  of  zeros  provided 
by  the  Jenkins  (1975)  program.  These  zeros  are  the  result  of 
factoring  the  frequency  response  provided  by  the  McClellan,  Parks 
and  Rabiner  (1973)  program. 

COMBO  first  generates  an  array  consisting  of  all  of  the  zeros 
(complex  conjugate  pairs  or  reals)  to  be  placed  in  Hj(z),  while  all 
others  are  assumed  to  be  placed  in  HgU).  The  program  insures  that 
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for  responses  which  contain  five  real  zeros,  a  sufficient  number  of 
combinations  occur.  Once  a  combination  has  been  specified,  program 
control  is  passed  to  the  zero  reconstitution  subroutine. 

The  reconstitution  subroutine  P0LYREC0N  is  responsible  for 
multiplying  the  appropriate  zeros  together  to  form  a  test  case 
Hj ( z )  and  H2(z) .  This  algorithm  uses  either  a  real  root  to  form 
a  linear  factor  or  a  complex  conjugate  pair  to  form  a  quadratic 
factor.  All  arithmetic  performed  is  real.  The  H^z)  and  H2(z) 
arrays  are  readily  converted  to  scaled  hj(n)  and  h2(n)  arrays. 

The  arrays  are  combined  as  described  in  the  FOM  procedure,  the  statis 
tics  are  performed  by  the  HSTAT  subroutine  and  an  FOM  is  assigned. 
Next,  the  FOM  is  compared  to  the  FOM  of  the  previous  design  and 
a  decision  is  made  as  to  which  is  best.  If  the  new  design  is  best, 
those  results  are  stored  and  the  program  proceeds  to  iterate  again. 
Upon  completion,  the  results  are  passed  back  to  the  calling  routine 
for  further  analysis. 


CHAPTER  VI 


COMPUTER  AIDED  DESIGN  APPLICATION 

The  Solid  State  Devices  Lab  group  at  the  University  of  Central 
Florida,  headed  by  Dr.  Donald  Malocha,  uses  a  computer  analysis 
system  called  SAWCAD,  developed  at  UCF,  to  design  SAW  filter 
devices.  Added  to  this  program  are  the  McClellan,  Parks  and 
Rabiner  (1973)  (Remez)  program,  the  Jenkins  (1975)  zero  location 
program  and  the  Optimizing  program  discussed  in  Chapter  V.  The 
following  is  a  brief  discussion  of  the  use  of  the  added  features 
to  SAWCAD.  The  complete  SAWCAD  package  is  not  discussed  here. 
Further  information  on  other  aspects  of  SAWCAD  can  be  obtained  by 
referring  to  Richie  (1983). 

A  listing  of  the  main  menu  is  shown  in  Figure  9.  A  filter 
design  is  initiated  by  selecting  the  (C)hebyshev  [REMEZ]  option  to 
the  main  menu.  The  program  proceeds  to  the  McClellan,  Parks  and 
Rabiner  (1973)  program  which  has  been  somewhat  modified.  The 
data  entry  routine  consists  of  an  interactive  session  between  the 
computer  and  the  designer.  During  this  session,  the  designer  is 
asked  to  specify  the  filter  type  (i.e.,  bandpass,  differentiator 
or  Hilbert  transformer),  the  number  of  distinct  frequency  bands  up 
to  fsamp-|e/2>  the  start  and  stop  frequencies  of  each  band,  the 
maximum  dB  level  in  each  band  and  the  maximum  allowable  ripple  in 
the  primary  passband  (multiple  passband  designs  are  possible). 


Figure  9.  SAWCAD  Main  Menu. 
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Once  frequency  range,  function  and  amplitude  information  is 
entered,  the  program  makes  an  initial  guess  of  the  required  filter 
order  and  initiates  a  design.  If  the  design  fails  to  converge  to 
the  required  specifications,  the  filter  order  is  increased  and  the 
process  is  repeated  until  convergence  is  obtained.  Once  a  valid 
design  is  found,  the  design  report  is  printed  and  control  is  passed 
back  to  the  SAWCAD  main  menu. 

At  this  point,  the  impulse  response  coefficients  exist  in 
memory  only.  The  SAWCAD  (W)rite  function  is  selected  and  a  disk 
file  containing  the  impulse  response  coefficients  may  be  written. 
This  step  is  highly  recomnended.  Control  passes  back  to  the  SAWCAD 
main  menu. 


Any  number  of  options  are  now  open  to  the  designer.  The 
obtained  response  could  be  used  inmediately  with  other  SAWCAD 
functions  if  desired  (i.e.,  FFT,  Graphics  analysis,  etc.).  The 
next  option  we  shall  concern  ourselves  with  here  is  the  (Z)ero 
Extraction  option.  Selecting  this  option  requires  no  further 
input,  since  the  program  calls  the  Jenkins  and  Traub  program  to 
factor  the  z-transform  of  the  impulse  response.  The  program 


returns  the  zeros  and  places  them  in  the  amp  and  phase  variables 
of  the  SAWCAD  program.  The  designer  must  be  aware  that  the  impulse 
response  is  no  longer  in  memory! 

Next,  the  optimum  split  for  low  order  designs  can  be  found 
using  the  (S)plit  Transducers  selection  on  the  SAWCAD  main  menu. 


No  further  user  input  is  required  since  the  program  iterates  until 
the  best  split  design  is  found.  Once  found,  the  program  prompts 
the  user  for  transducer  1  and  2  file  names  under  which  to  save  the 
impulse  response  coefficients. 

Generation  of  the  design  is  now  complete.  It  may  be  checked 
by  multiplying  the  transducer  1  and  2  data  files  together  and 
comparing  the  product  with  the  original  response  generated  by  the 
Remez  method.  These  comparisons  may  be  done  graphically  using  the 
powerful  graphics  and  FFT  facilities  of  the  SAWCAD  environment. 


CHAPTER  VII 


RESULTS 

In  an  attempt  to  determine  a  general  algorithm  applicable  to 
filters  of  any  order  or  type,  eight  different  low-order  test  filters 
were  designed.  Four  of  the  filters  were  low  pass  designs  and  four 
were  high  pass.  The  input  design  specs  appear  in  Table  3. 

The  designs  were  made  using  the  Remez  technique,  factored  by 
the  Jenkins  (1975)  program,  and  were  subjected  to  the  optimizing 
program  to  test  all  possible  split  designs.  The  composition  of  each 
transducer  was  recorded  to  reflect  the  number  of  stopband  and  passband 
zeros,  and  whether  these  zeros  were  real  or  complex.  These  results 
are  tabulated  in  Table  4. 

Careful  study  of  the  results  does  not  indicate  any  clearly 
emerging  pattern.  In  three  out  of  four  of  the  cases  with  very 
narrow  passbands  (LP1,  HP1  and  HP4),  the  passband  zeros  all  appeared 
on  one  transducer  accompanied  by  several  stopband  zeros.  However, 

LP4  passband  zeros  were  split  between  the  two  transducers.  The 
other  cases  did  not  produce  results  which  might  indicate  a 
predictable  pattern. 

Another  test  was  devised  to  test  and  rate  a  conventional 
no-split  design,  a  strictly  passband-stopband  split,  the  "alternating 
zero"  algorithm  employed  by  other  design  groups  (Morimoto  et  al . 


stopband,  P  =  passband 
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1980  and  Ruppel  et  al.  1984),  and  the  split  algorithm  presented 
here.  The  HP1  filter  spec  was  arbitrarily  chosen  as  the  subject  of 
the  test.  Figure  10  shows  the  overall  frequency  response  of  the 
HP1  filter. 

The  conventional  no-split  design  yielded  a  figure  of  merit 
(FOM)  of  about  2.5214,  whereas  the  passband-stopband  split  (i.e., 
all  passband  zeros  on  one  transducer  and  all  stopband  zeros  on  the 
other)  yielded  an  FOM  of  about  3.588.  Figures  11  and  12  show  the 
frequency  response  of  each  transducer.  The  passband  transducer 
shows  nearly  unity  gain  at  0  with  a  gradual  rolloff  in  the  region 
of  f  ,  a  very  difficult  response  to  implement  on  a  SAW  device. 

The  "alternating  zero"  approach  yielded  a  FOM  of  about  6.5756. 
This  shows  an  improvement  in  the  quality  of  the  design  as  compared 
to  the  passband-stopband  split  method.  Figures  13  and  14  show  the 
frequency  response  of  each  transducer.  The  zeros  are  clearly 
orthogonal.  Figures  15  and  16  show  the  impulse  response  represen¬ 
tations.  This  design  does  appear  to  have  merit  in  the  case  where 
a  fifty-fifty  split  is  highly  desired. 

The  technique  presented  here  produced  a  design  with  a  FOM  of 
about  9.753.  Figures  17  and  18  show  the  frequency  response  of  each 
transducer  and  figures  19  and  20  show  the  corresponding  impulse 
response  plots. 

The  last  two  designs  do  not  appear  to  pose  fabrication  problems 
which  might  plague  the  passband-stopband  split  design  case.  As  a 


check  of  the  split,  the  frequency  responses  of  the  last  design  were 
multiplied  together  via  SAWCAD  and  replotted.  That  plot  is  exactly 
the  same  as  the  original  frequency  response  in  Figure  10,  as 
expected. 


Figure  10.  Overall  Frequency  Response. 


Figure  11.  Transducer  1  Frequency  Response  (Passband-Stopband  Design). 


Transducer  2  Frequency  Response  (Passband-Stopband  Design) 


Transducer  1  Frequency  Response  (Alternating  Ze 


Figure  14.  Transducer  2  Frequency  Response  (Alternating  Zero  Design). 


Transducer  1  Impulse  Response  (Alternating  Zero  Design) 


Transducer  2  Impulse  Response  (Alternating  Zero  Design) 


Figure  17.  Transducer  1  Frequency  Response  (Fully  Optimized  Design). 


Figure  19.  Transducer  1  Impulse  Response  (Fully  Optimized  Design) 


CHAPTER  VIII 


CONCLUSIONS 

The  technique  presented  in  this  thesis  of  optimally  splitting 
the  zeros  of  the  transfer  function  shows  considerable  immediate 
promise  for  low  order  (N  less  than  40)  filter  designs.  An  algorithm 
permitting  the  split  design  of  high  order  filters  has  not  become 
readily  apparent,  although  the  presented  "all -combi nations"  technique 
may  still  be  used  if  limits  are  placed  upon  the  transducer  sizes. 

For  example,  specifying  a  fifty-fifty  split  design  significantly 
reduces  the  number  of  combinations  which  must  be  tested.  Future 
efforts  in  this  area  must  focus  upon  a  more  efficient  search  criteria 
or  be  content  to  accept  less  than  optimal  results,  as  in  the 
"alternating  zero"  approach.  The  concepts  presented  here  may  be 
extended  to  generate  multi -transducer  splits,  with  primary  utility 
in  large  ordered  digital  filter  implementations. 

In  spite  of  the  limitations  of  the  technique  used  here,  it  is 
obvious  that  a  "best  split"  design  of  an  FIR  filter  transfer 
function  exists  based  upon  the  average  size,  variance  and  range  of 
the  resulting  impulse  response  coefficients.  The  "alternating  zero" 
split  approach  (Morimoto  et  al.  1980  and  Ruppel  et  al .  1984),  the 
only  other  method  published  to  date,  exhibited  a  lower  FOM  when 
compared  to  the  split  found  by  testing  all  possible  combinations, 
based  upon  the  stated  criteria.  Careful  study  of  figures  15,  16,  19 


71 


1 033C3R3? 


<5 


72 


and  20  (shown  in  Chapter  VII)  demonstrates  that  the  presented 
technique  does,  indeed,  produce  larger  tap  sizes  for  a  given  transfer 
function  than  the  "alternating  zero"  approach.  Therefore,  for 
critical,  low-ordered  design  cases,  the  technique  presented  here 
will  provide  best  results. 

The  Jenkins  and  Traub  factoring  algorithm  is  a  very  accurate 
means  of  obtaining  the  zeros  of  polynomials  within  an  order  of  one- 
hundred  or  so.  However,  high  ordered  SAW  filter  designs  will  require 
that  the  zeros  of  polynomials  of  a  thousand  order,  or  more,  will  need 
to  be  extracted.  Based  upon  observation  of  the  programs  and 
polynomial  characteristics  encountered  in  this  thesis,  several 
alternative  approaches  to  locating  the  zeros  seem  feasible. 

One  technique  that  merits  further  investigation  is  the  use  of 
the  Fourier  transform  to  expose  the  stopband  zero  locations.  This 
can  be  done  by  noting  where  the  stopband  nulls  occur  and  reference 
these  points  to  the  corresponding  angles  about  the  z-plane  unit 
circle.  The  cosines  and  sines  of  these  angles  are  the  real  and 
imaginary  components  of  the  stopband  zeros  which,  of  course,  have 
corresponding  complex  conjugates  in  the  lower  half  of  the  unit  circle. 
Once  all  of  the  stopband  zeros  have  been  found,  the  z-polynomial  may 
be  deflated  in  one  step,  yielding  a  polynomial  consisting  of  only 
the  passband  zeros.  Next,  the  Jenkins  and  Traub  algorithm  may  be 
applied  to  the  greatly  reduced  polynomial  to  locate  the  passband 
zeros.  The  Fourier  transform  technique  could  be  applied  to  any  FIR 
obtained  by  any  design  technique. 


Another  zero  location  technique  which  may  be  explored  is  unique 
to  the  Parks  and  McClellan  program.  One  of  the  outputs  of  the 
program  is  a  listing  of  all  of  the  Chebyshev  polynomial  extremal 
frequencies.  Since  the  limits  of  each  stopband  are  known  (specified 
apriori),  and  since  the  extremal  frequencies  are  fairly  evenly 
spaced  in  the  stopband,  an  interpolation  between  adjacent  extremal 
frequencies  will  provide  a  good  approximation  (or,  at  least  a  good 
initial  guess)  to  a  function  zero  in  the  stopband.  If  the 
interpolation  method  is  deemed  sufficient,  then  all  of  the  zeros 
found  may  be  used  to  reduce  the  response  down  to  a  polynomial 
containing  only  the  passband  zeros,  to  which  the  Jenkins  and  Traub 
algorithm  may  be  applied.  A  better  approximation  of  the  stopband 
zeros  would  be  obtained  by  applying  the  interpolation  approximation 
as  an  initial  guess  to  Bairstow's  technique  to  obtain  the  best 
estimate  of  the  zero.  Again,  the  passband  zeros  could  be  found 
using  the  Jenkins  and  Traub  method. 

In  conclusion,  an  optimal  split  of  zeros  does  exist,  based  upon 
the  stated  criteria.  The  major  limiting  factor  of  the  presented 
technique  is  the  requirement  that  all  possible  split  combinations 
must  be  devised  and  rated  with  respect  to  each  other.  This  is  a 
very  time-consuming  process  and  future  efforts  may  lead  to  a  more 
efficient  means  of  finding  the  optimal  combination.  The  use  of  more 
powerful  computing  machines  may  make  the  split-design  of  somewhat 
higher  order  filters  practical,  but  a  reasonable  upper  limit  is 


rapidly  approached  with  each  increase  in  filter  order.  Perhaps  the 
best  tradeoff  between  this  method  and  the  "alternating  zero"  approach 
is  to  obtain  the  best  fifty-fifty  split  by  trying  all  possible 
combinations  for  this  case  only  (i.e.,  N  zeros  taken  N/2  at  a  time). 
This  approach  will,  at  the  very  least,  provide  as  good  a  design  as 
the  "alternating  zero"  approach,  without  the  serious  high  order 
computational  drawbacks  of  the  method  presented  here.  Additionally, 
some  design  considerations  may  favor  a  fifty-fifty  split,  as  in  the 
case  of  two  digital  signal  processors  being  required  to  evenly  share 
the  processing  tasks.  Nevertheless,  the  purpose  of  this  thesis  was 
to  prove  the  existence  of  one  such  optimal  combination,  not  the 
optimal  means  of  obtaining  it.  With  a  sigh  of  relief,  and  a  hint 
of  moderate  surprise,  such  proof  has  been  presented. 
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PROGRAM  FOR  THE  DESIGN  OF  LINEAR  PHASE  FINITE  IMPULSE 
RESPONSE  (FIR)  FILTERS  USING  THE  REMEZ  EXCHANGE 
ALGORITHM 

C  JIM  MCCLELLAN,  RICE  UNIVERSITY,  APRIL  13,  1973 
C 

C  MODIFIED  BY  KEITH  V.  LINDSAY,  UNIVERSITY  OF 
C  CENTRAL  FLORIDA,  1  MARCH  86,  FOR  USE  WITH  UCF'S 
C  SflWCAD  PROGRAM. 

C 

c  THREE  TYPES  OP  FILTERS  ARE  INCLUDED— BANDPASS  FILTERS 
C  DIFFERENTIATORS,  AND  HILBERT  TRANSFORM  FILTERS 
C 

C  THE  INPUT  DATA  CONSISTS  OF  5  SECTIONS 
C 

C  SECTION  1— FILTER  LENGTH,  TYPE  OF  FILTER,  1-MJLTIPLE 
C  PASSBANiy STOPBAND,  2-DIFFERENTIATOR,  3 -HILBERT  TRANSFORM 
C  FILTER,  NUECER  OF  BANDS,  CARD  PUNCH  DESIRED,  AND  GRID 
C  DENSITY 
C 

C  SECTION  2— BANDEDGES,  LOWER  AND  UPPER  EDGES  FOR  EACH 
BAND 

C  WITH  A  MAXMJM  OF  10  BANDS. 

C 

C  SECTION  3— DESIRED  FUNCTION  (OR  DESIRED  SLOPE  IF  A 
C  DIFFERENTIATOR  )  FOR  EACH  BAND. 

C 

C  SECTION  4— WEIGHT  FUNCTION  IN  EACH  BAND.  FOR  A 
C  DIFFERENTIATOR,  THE  WEIGHT  FUNCTION  IS  INVERSELY 
C  PROPORTIONAL  TO  F 
C 

OOMMON/FILE/  AMP  (40 96)  ,PHASE(4096)  ,NFFT,  ITYPE 

COMMON/ DAT/  FO,  TFLO,  TFHI,  NUM 

COMMON 

PI2 ,  AD,  DEV,  X,  Y,  GRID,  EES,  WT,  ALPHA,  IEXT,  NFCNS,  N3RID 

DIMENSION  IEXT(514)  ,AD(514)  ,ALPHA(514)  ,X(514)  ,Y(514) 
DIMENSION  H(514) 

DIMENSION  EES  (8224)  ,GRID(8224)  ,WT(8224) 

DIMENSION  EDGE (20)  ,FX(10)  ,WTX(10)  ,DEVIAT(10) 

DOUBLE  PRECISION  PD2,PI 
DOUBLE  PRECIS  ION  AD,DEV,X,Y 
DOUBLE  PRECISION  HB 
LOGICAL  FAIL,  MOREN 
INTEGER  PB,  SB 
PI2^  .283185307179586 
PI-3 .1415926535897  S3 
C 

C  THE  PROGRAM  IS  SET  UP  FOR  A  MAXIMUM  LENGTH  OF  1024,  BUT 
C  THIS  UPPER  LIMIT  CAN  BE  CHANGED  BY  REDIMENSICNING  THE 
C  ARRAYS  IEXT,  AD,  ALPHA,  X,  Y,  H  TO  BE  NEMAX/2+2. 

C  THE  ARRAYS  EES,  GRID,  AND  WT  MUST  BE  DIMENSIONED 


C  16(NEMAX/2+2) . 

C 

NFMAX=1024 
100  CONTINUE 
JTXI&O 
C 

C  PROGRAM  INHIT  SECTION 
C 

HUNT  *, '  DIGITAL  FILTER  DESIGN  (FIR)  VIA  THE 
1  REMEZ  EXCHANGE  ALGORITHM.  1 
PRINT  '  ' 

HUNT  *,  ’ENTER  TYPE  OF  FILTER: 1 

HUNT*,'  (1) -MULTIPLE  PASEBAND/STDPBAND' 

HUNT  *, '  (2)  -DIFFERENTIATOR' 

PRINT  *, '  (3)  -HILBERT  TRANSFORM  FILTER' 

PRINT  ' 

READ  *,  JTYPE 
PRINT  ' 

PRINT  *,  'ENTER  THE  NUMBER  OF  BANDS  ' 

READ  *,IBANDS 
PRINT  • 

C  PRINT  *,  'OUTPUT  THE  IMPULSE  RESPONSE  (1=YES,  0=NC 
C  READ  *,JPUNCH 
JPUNCB=1 

C  PRINT  *,  'ENTER  THE  GRID  DENSITY' 

C  READ  *fIGRID 
LGR3Df16 

IF(FBANDS.LE.O)  EBANE6=1 
C 

C  GRID  DENSITY  IS  ASSUMED  TO  BE  16  UNLESS  SPECIFIED 
OTHEIWISE. 

C 

IF  (LGRID.  LE.  0)  IGRDXL6 
DO  888  JKL,fBANDS 
PRINT  V  ' 

PRINT  *,  'BAND  ' ,  J, ' : ' 

PRINT  *, '  POWER  EDGE: ' 

READ  *f  EDGE(2M-1) 

PRINT  V  ' 

PRINT  *, '  UPPER  EDGE: ' 

READ  *,EDGE(2\T) 

888  CONTINUE 
IF(JTYPE.EQ.2)  GO  TO  890 
PRINT  V  ' 

PRINT  *, '  ENTER  THE  DESIRED  FUNCTION  OF  EACH  BAND' 
PRINT  *, '  (O-NOEASS,  1-PASffiAND)  • 

DO  889  i>l,FBAND6 
PRINT  ' 

PRINT  *,  ’  BAND  ' ,  J, ' : ' 

READ  *,FX(J) 

889  CONTINUE 
GO  TO  893 

890  DO  892  *>1,IBAND6 
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HUNT  *, 'ENTER  THE  SLOPE  OF  BAND  1  ,J 
READ  *,FX(J) 

892  CONTINUE 

893  CONTINUE 
HUNT  V  ' 

HUNT  *,  'ENTER  THE  MAX  cB  RIPPLE  IN  THE  EASffiAND' 

READ  *,DP 

DP°10.  ED**  (DP/20)  -l.OBO 
DS-l.OED 

DO  898  J»1,*BAND6 
HUNT  *,'  ' 

HUNT  *,  'ENTER  THE  MAXIMUM  dB  LEVEL  IN  BAND  ' , J 
READ  *,WTX(J) 

IF  (WTX(J)  .B3.O.0E0)  THEN 
WTX(J)  =1.0  ED 

GO  TD  898 

END  IF 

C  WTX(J)  =AINT(10.0ED**(  (20.0ED*ALOG10(DP)  - 
WTX(J))/20.0ED)  +1) 

WTX(J)  »(10.0EO**  ( (20.0ED*ALCG10(DP)  -WTX(J) )  /20.0ED) ) 
IF  (WTX(J)  .LT.l.OED)  WTX(J) «1.0ED 
HUNT  *,  'STOEBAND  WEIGHTING  =',WTX(J) 

2477  IF  (DS.GT.EP/WTX(J))  THEN 

DS*=DP/WTX(J) 

IF  (J.EQ.l)  THEN 

EELTAF«EDGE(3)  -EDGE  (2) 

GO  TO  898 

END  IF 

IF  (J.  LT.  JBANDS)  THEN 

IF  (DP/ViTX(J-l)  .EQ.DP)  THEN 

CELTAF»EEGE(J*2-I)  -EEGE(J*2-2) 
GO  TO  898 

END  IF 

DELTAF«EDGE(J*2+1)  -EDGECJ*2) 

GO  TO  898 

END  IF 

IF  (J.EQ.EBANDS)  THEN 

EELTAF=ELGE(J*2-1)  -EDGE(J*2-2) 

END  IF 

END  IF 
898  CONTINUE 
C 

C  TAKE  A  GUESS  AT  AN  INITIAL  VALUE  FOR  NFILT 
C  BASED  UPON  VAIDaNATHAN'S  FORMULATION 
C 

NFILTV2NT((-10.0ED*AL0G10(DP*DS)  -  13.0ED) /(14.6ED  * 
DELTAF) ) 

2126  HUNT  *,  'WORKING  CN  FILTER  OF  ORDER  '  ,NFILT-1 

IF (NFILT.GT.NEMAX.CR.NETLT.LT.  3)  CALL  ERROR 
IF(JHPE.EQ.O)  CALL  ERROR 
NB>1 
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IF  (JTYPE.  EQ.l)  NB3-0 
NCCD=NFILT/2 
NCCONFILT-2  *NCCO 
NFCN9=NFILT/2 

IF(NCDD.HM.AND.NB3.H3.0)  NPCN9=NPCNSH 
C 

C  SET  UP  THE  CENSE  GRID.  THE  NUMBER  OF  POINTS  IN  THE  GRID 
C  IS  (FILTER  LENGTH  +  1)  *GRID  DENSITY  /  2 
C 

GRID(l)  =EEGE(1) 

DELF=IGRID*NPCNS 
EELF=0 .5/DELF 
IF(NBG.EQ.O)  GO  TO  135 
IF  (EDGE  (1)  .LT.DEILF)  GRID  (1)  =>€ELF 
135  CONTINUE 
J=1 
L=1 

IBAND-1 

140  FUI^EDGE(LH) 

145  TEMP»GRID(J) 

C 

C  CALCULATE  THE  DESIRED  MAGNITUDE  RESPONSE  AND  THE  WEIGHT 
C  FUNCTION  CN  THE  GRID 
C 

EES  ( J)  -EFF  (TEMP,  FX,  WTX,  IBAND,  JTYPE) 

WT(J)  =WATE  (TEMP,  FX,  WTX,  IBAND,  JTYPE) 

CKH1 

GRID(J)  =TEMP4-DELF 
IF  (GRID  (J)  .GT.FUP)  GO  TO  150 
GO  TO  145 
150  GRID(J-l)  -FUP 

DES  (J-l)  -EFF  (FUP,  FX,  WTX,  IBAND,  JTYPE) 

WT(J-l)  -WKTE  (FUP,  FX,  WTX,  IBAND,  JTYPE) 

IfiANIMBANDfl 

>Lf2 

IF  (IBAND.  GT.  WANDS)  GO  TO  160 
GRID  (J) -EDGE  (L) 

GO  TO  140 
160  IGR33XJ-1 

IF  (NB3.  NE.  NOCD)  GO  TO  165 
IF(GRID(!GRID)  .GT.  (0.5-DELF) )  N3RID-N3RID-1 
165  CONTINUE 
C 

C  SETT  UP  A  NEW  APERCKEMATICN  HKBLEM  WHICH  IS  H2DIVALENT 
C  TO  THE  ORIGINAL  IRCBLEM 
C 

IF(NHj)  170,170,180 
170  IF(N0CD.H3.1)  GO  TO  200 
DO  175  i>1,NjRID 
CHANG  B-DCDS  (PI  *GRID(J) ) 

DES  ( J)  «OES  ( J)  /CHANGE 
175  WT(J)^iT(J)*<HANGE 
GO  TO  200 
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180  IF(NCCD.BQ.l)  GO  TO  190 
DO  185  i>l,N3RID 
(EAHGE=DSIN(PI<GRID(J) ) 

DES  (J)  -OES  (J)  /CHANGE 
185  WT(J)  -WT(J)  *CHANGE 
GO  TO  200 

190  DO  195  >>1,JGRID 

CHAN3E>=DSIN(PI2*GRID(J) ) 

DES(J) -DES(J) /CHANGE 
195  WT(J)  ^T(J)  ^CHANGE 
C 

C  INITIAL  GUESS  FOR  THE  EXTREMAL  ERR2UENCIES — EQUALLY 
C  SPACED  ALONG  THE  GRID 
C 

200  TEJO^FLQAT  (EGRID-1)  /FLOAT (NPCNS) 

DO  210  J“1,NFCNS 
210  IEXT(J)  =(J-1)  *TEWPfl 
IEXT  (NPCNSf  1)  “NS  RID 
NM1-NPCN&-1 
JE«NPCN9fl 
C 

C  GALL  THE  REHEZ  EXCHANGE  ALGORITHM  TO  DO  THE 
APPROXIMATION 
C  PROBLEM. 

C 

CALL  REMEZ  (EDGE,  EBANDS,  MOREN) 

IF  (M3REN.  EQ . .  TKJE. )  THEN 
NFILTVNFILTH 
GO  TO  2126 

END  IF 

IF  <DEV/WTX(PB)  .GT.DP  .OR.  EEV/WTX(ffi)  .GT.DS)  THEN 
NFTLTVNFILTH 

C  PRINT  *, 'DEVIATION  PB  =  ' ,DEV/WIX(PB) , 'IB- '  ,IB 

C  PRINT  V DEVIATION  SB  ,EEV/WEC<ffl) 

GO  TO  2126 

END  IF 
C 

C  CALCULATE  THE  IMPULSE  RESPONSE. 

C 

IF  (NEE)  300,300,320 
300  IF(NQCD.EQ.O)  GO  TO  310 
DO  305  i>l,N£L 
305  H( J)  .5*ALPHA(EC-J) 

H(NFCNS)-ALPHAd) 

GO  TO  350 

310  H(l)  >0 .25*ALPHA(NFCNS) 

DO  315  Cb2,  Ml 

315  H(J) -0.25* (ALPHA (IE-J) 4ALPHA(NFCNSf2-J) ) 
H(NFCNS)*0.5*ALPBA(1)  -K).25*ALPHA(2) 

GO  TO  350 

320  IF(NOCD.BQ.O)  GO  TO  330 
H(1M>  .25*ALPHA(NFCNS) 

H(2)^)  .25*ALPHA«ML) 


DO  325  >3, WO. 

325  H< J)  ^)  .25*  (ALPHA  (EE -J)  -ALPHA  (NFCNS*-3-J) ) 

H(NPCNS)  -0 .5*ALPHA<1)  -0 .25*ALPHA(3) 

H(1E)«0.0 
GO  TO  350 

330  H(1)^).25*ALPHA(NFCNS) 

DO  335  J>2,!KL 

335  H(J)  =0.25*(ALPHA(IE,-J) -ALPHA  (NFCNSf2-J) ) 

H(NPCNS)  -0 .5*ALPHA(1)  -0 .25*ALPHA(2) 

C 

C  SET  UP  IMPULSE  RESPONSE/POLYNOMIAL  ARRAY  ** 

C  ARRAY  IN  VARIABLE  AMP  ** 

C 

350  DO  342  I»1,NFCNS 
AMP(I)  =H(I)  /H(NFCNS) 

IP(NB3.«J.0)  AMP(NFILT-I+1)«H(I)/H(NFCNS) 

IF(NHG.EQ.l)  AMP(NFILTfl-I)— H(I)/H(NFCNS) 

342  CONTINUE 

IFCNB3.EQ.1  .AND.  N0DD.EQ.1)  AMP(fE)O.DO 
C 

C  ADD  SWCAD  PARAMETERS  NORMALIZED  TO  1  MHZ 
C  2  Fo  SAMPLING 
C 

NUfrNFILT 

NFETVNFILT 

ITYPE>-1 

FCN0.5 

TELO-(NFILT-l)/2 
TSHI«  (NFILT-1)  /2 
C 

C  PROGRAM  OUTPUT  SECTION. 

C 

C 

PRINT  360 

360  FORMAT (//70(1H*)//25X,  'FINITE  IMPULSE  RESPONSE  (FIR)  V 

1  25X,  'LINEAR  PRASE  DIGITAL  FILTER  DESIGN'/ 

2  25X,'REMEZ  EXCHANGE  ALGORITHM'/) 

IF(JTHE.EQ.l)  PRINT  365 

365  FORMAT (25X,  'BANDPASS  FILTER'/) 

IF(JTOPE.E3.2)  PRINT  370 
370  FORMAT (25X,  'DIFFERENTIATOR'/) 

IF(JTZPE.EQ.3)  PRINT  375 
375  FORMAT (25X,  '  HILBERT  TRANSFORMER'/) 

PRINT  378rNFILT 

378  FORMAT (20X,  'FILTER  LENGTH  -  '  ,13/) 

PRINT  380 

380  FQRMAT(20X# '  *****  IMPULSE  RESPONSE  ***♦*') 

DO  381  J«l,NFOJS 

K-NFILTfl-J 

IF(NB3.E3.0)  PRINT  382, J,H(J),K 
IF(NHG.E3.1)  PRINT  383, J, B(J)  ,K 

381  CONTINUE 

382  F0RMAT(20X,'H(M3,')  -  ',EL5.8,'  -H(M4,')') 


83 

383  FQRMAT(20X,  *H( ' *13* ')  -  *,E15.8,  '«  -H(  ',14,')  •) 
IF(NEG.EQ.l.AND.NODD.EQ.l)  PRINT  384, IE 

384  FORMAT (2GX,'H<',  13,')  -O.O') 

DO  450  R=l,fBANIS,4 
KUP=K+3 

IF  (KUP.  GT.  N3AND6)  KUE*JBAND6 
PRINT  385,<J,J»K,KUP) 

385  FQRMAT(/24X,4(  'BAND'  ,I3,8X)) 

PRINT  390,  (EDGE(2*J-1)  ,J-K,KUP) 

390  PQRMAT(2X,  'LOWER  BAND  EDGE1 ,5F15.9) 

PRINT  395,  (EDGE(2^J)  ,v>K,KUP) 

395  FORMAT (2X,  'UPPER  BAND  EDGE*  ,5F15.9) 

IF (JTYPE.NE.2)  PRINT  400,  (FX(J)  „>K,HUP) 

400  F0RMAT(2X,  'DESIRED  VALUE1 ,2X,5F15.9) 

IF(JKPE.EQ.2)  PRINT  405,  (FX(J) , J*K,KDP) 

405  FQRMAT(2X,  'DESIRED  SLOPE' ,2X,5F15.9) 

PRINT  410,(WTX(J)  ,J*=K,KDP) 

410  EORMAT(2X,  'WEIGHTING1 ,6X,5F15.9) 

DO  420  J>K,RUP 
420  DEVIAT(J)  -DEV/WTX(J) 

PRINT  425,  (DEVIAT(J)  ,>K,KUP) 

425  FQRMAT(2X, 'DEVIATION' ,6X,5F15.9) 

IF  (JTYPE.  NE.l)  GO  TO  450 
DO  430  J=K,KUP 

430  DEVIAT(J)  =20.0*ALOG10(DEVIAT(J) ) 

PRINT  435,  (DEVIAT  ( J)  ,J-K,KUP) 

435  PQRMAT(2X,  'DEVIATION  IN  DB*  ,5F15.9) 

450  CONTINUE 

PRINT  455,  (GRID(IEXKJ) )  ,J=L,EE) 

455  FORMAT(/2X, '  EXTREMAL  EREQUQTCIES'/(2X,5F12.7) ) 
PRINT  460 

460  FQRMAT(/1X,70(1H*)/1H1) 

C 

RETURN 

END 

C 

EUNCTIDN  EFF  (TEMP,  FX,  WTX,  IB  AND,  JTYPE) 

c 

c  function  to  calculate  the  desired  magnitude  response 
c  as  a  function  of  frequency, 
c 

DIMENSION  FX(5)  ,WTX(5) 

IF(JT5fPE.E0.2)  GO  TO  1 
EFFHFXdBAND) 

RETURN 

1  EFF*FX(IBAND)  *TEMP 
RETURN 
END 


c 

c 

c 


FUNCTION  WATE  (TH1P,  FX,  WJX,  IBAND,  JTYPE) 


non 
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c  function  to  calculate  the  weight  function  as  a  function 
c  of  frequency, 
c 

DIMENSION  FX(5)  ,WTX(5) 

IF(JTYPE.  EQ.2)  GO  ID  1 
WATE>=WTX(I£AND) 

RETURN 

1  IF (FX(IBAND)  .LT. 0.0001)  GO  TO  2 
WATE>=WTX  (IBM©)  /TEMP 

RETURN 

2  WATE*WTX(I£AND) 

RETURN 

END 

C 

C 

C 

SUBROUTINE  ERROR 
HUNT  1 

1  FORMAT  ( '  **********  ERROR  IN  INPUT  DATA  ***********) 
STOP 
END 

c 

c 

c 

SUBROUTINE  REMEZ  (EDGE,  1BANE6,  M3REN) 

C 

c  this  subroutine  implenents  the  remez  exchange  algorithm 
c  for  the  weighted  chetycher  approximation  of  a  continuous 
c  function  with  a  sum  of  ao sines,  inputs  to  the 
subroutine 

c  are  a  dense  grid  which  replaces  the  frequency  axis,  the 
c  desired  function  on  this  grid,  the  weight  function  on 
this 

c  grid,  the  nunber  of  cosines,  and  an  initial  guess  of  the 
c  extranal  frequencies,  the  pcogran  minimizes  the 
chefcychev 

c  error  by  determining  the  best  location  of  the  extrenal 
c  frequencies  (points  of  maxi, min  error)  and  then 
calculates 

c  the  coefficients  of  the  best  approximation. 

OMEN 

PI2  ,AD,EEV,  X,  Y,  GRID,  EES,WT,  ALPHA,  TEXT,  NFCNS,  N3RID 
DIMENSION  EDGE  (20) 

DIMENSION  XEXT(514)  ,AD(514)  ,ALEHA(514)  ,X(514) ,  Y(514) 
DIMENSION  EES (8224) , GRID (82 2 4)  ,WT(8224) 

DIMENSION  A(514)  ,P(513)  ,Q(513) 

DOUBLE  PRECISION  PI2,ENUM,ECEN,DTEMP,A,P,Q 
DOUBLE  DECISION  AD,EEV,X,Y 
LOGICAL  MDREN 

IBE  PROGRAM  ALLCWS  A  MAXIMUM  NOTBER  OF  ITERATIONS  OF  25 


M3RH*. FALSE, 


nSMAX-25 
DEVI>-1.0 
tE»NEOJSfl 
ICZ*NFCNS+2 
NnEfcO 
100  CONTINUE 

IEXT  (IB  2)  =N3RU>1 
NITEF*=NITERfl 

IF  ( NITER*  GT.  ITOMAX)  GO  TO  400 
DO  110  .>1,1E 
DOTEMP*GRID(IEXr(J) ) 

DTEKWXDS(ETEMP*PI2) 

110  X(J)-D!TEMP 

JEIb(NPCN9-l)/15+l 
DO  120 

120  AD(J)>=D(J,IB,JET) 

mufro.o 

EDE»«0.0 

K-l 

DO  130  J-1,NZ 
I^IEXr(J) 

DTHB*AD(J)  *DES(L) 

ENUFkENUHt'DTEMP 
nTE30^K*AD(J)  /WT(L) 

DDEN-DOEN+DTEMP 
130  K— K 

DEV-ENUM/DCEN 

ND*1 

IF  (DEV.  GT.  0.0)  NO— 1 

DEV—ND«DEV 

K-ND 

DO  140 

I>IEXr(J) 

DTEMM(*DEV/WT(L) 
y(j)"DES(L)+nrmp 
140  K»  K 

IF(DEV.GE.EEVL)  GO  TO  150 
MORJN-.TRDE. 

RETUR1 
150  DEVD-DEV 
JCHMGB-0 
Kl-IEXr(l) 

KNZ-IEXr(tE) 

KLCW-0 
NUT— ND 
J-l 
C 

C  SEARCH  FOR  THE  EHHEMAL  ERH2UINCIES  OF  TEE  BEST 
C  AFERCDOMATION 
C 

200  IF(J.Ba.ICZ)  BB-OOMP 
IF(J.GE.JBZ)  GO  TO  300 
KDMEXT(J+1) 


]>iExr(j)+i 

NOT^-NOT 

IF(J.EQ.2)  Yl=ODMP 
CDKE^EEV 

IF(L.GE.KUP)  GO  TO  220 
EWMIEE(L,IE) 

ERR-  (ERR-DES(L) )  %T(L) 
nraiwcT*ERR-coMP 
IF(nTEMP.LE.O.O)  GO  TO  220 
OOI!P»NOT*ERR 
210  J>Lfl 

IF(L-GE.RQP)  GO  TO  215 
ERR-GEE(L,JE) 
ERI*(ERR-DES(L))%T(D 
ETEMMBT*ERR-CDMP 
IF(DTEMP.LE.0.0)  GO  TO  215 
03MP«N0T*ERR 
GO  TO  210 
215  IEXT(J)  -Ir-1 
J=Jfl 
KLCW-L-1 
JCHJCBsJCHIGEfl 
GO  TO  200 
220  Ij-Ir-1 
225  D-I/-1 

IF  (L.  LE.  KLCW)  GO  TO  250 
&%R-GEE(L,lO 
ERR‘(ERR-EES(L})*WT(L) 
E(PEMW17P*ERR-CDMP 
IF(ETEMP.GT.O.O)  GO  TO  230 
IF  ( JGQEE.  LE.  0)  GO  TO  225 
GO  TO  260 
230  OOMP-NOT*ERR 
235  I>Ir-l 

IF(L.LE.KLCW)  GO  TO  240 
JsA*Gl!iE(LrI&} 
ERR-(ERR-DES(L))«WT(L) 
ETEMWHT*ERR-<X)MP 
IF  (ETTEMP.  LE.  0 .0)  GO  TO  240 
OQMP*M7T*ERR 
GO  TO  235 
240  KLCW-IEXr(J) 

IEXT(J)*Lfl 

J«Jfl 

JQUGEKJGE&GBfl 
GO  TO  200 
250  I*IEXr(J)+l 

IF(JCHICE.GT.O)  GO  TO  215 
255  I^Lfl 

IF(L.GE.KDP)  GO  TO  260 
ERIK3EE(L,ie) 

ERR*  (ERR- EES  (L) )  *WT(L) 
nTE^NUT*ERR-OOMP 


8; 


IP(nrEMP.LE.O.O)  GO  TO  255 
(30MP*NUT*ERR 
GO  TO  210 
260  KLCW=IEXT(J) 

J*J*1 
GO  TO  200 

300  IP(J.GT.NZZ)  GO  TO  320 

IF(Kl.GT.IEXra))  K1=IEXT(1) 
IP(KNZ.LT.IEXraC))  KNZ*IEXP(NZ) 

Nun-njr 

Ntn^-NO 

lr=0 

KUP-K1 

OOMMJZ*  (1.00001) 

LXK«1 
310  Irfrf-1 

IP(L.GE.KUP)  GO  TO  315 
ERIMjEE(LrlC) 
ERR=(ERR-DES(L))*WT(L) 
ETEME^NOT^ERR-OOHP 
IF(DTEMP.LE.O.O)  GO  TO  310 
OOM**11JT*ERR 
J«MZ 
GO  TO  210 
315  UJCK=6 
GO  TO  325 

320  IP  (LOCK.  GT. 9)  GO  TO  350 
IP(OQHP.GT.yi)  Y1KDHP 
Kl-IEXr(NSZ) 

325  MGRXDtl 
KLCW-KNZ 
NCnV-MTEL 
OOMI^Yl* (1.00001) 

330  Iflrl 

IP(L.LE.KLCW)  GO  TO  340 

HIEM5EE(L,)C) 

ERI*(ERR-DES(L))*fT(L) 

OTE3^N0T*ERR-O)MP 

IP  (OTEMP.LE.0 .0)  GO  TO  330 

j*iez 

OOMP»NOT*ERR 
LOO(«LOCK+10 
GO  TO  235 

340  IF  (LOCK.  DQ.  6)  GO  TO  370 
DO  345  <>1,NFCNS 
345  IEXT(ICZ-J)  -mr(ic-j) 

IEXT(1)-K1 
GO  TO  100 
350  IQ^IEXrUEZ) 

DO  360  J-1,NFCNS 
360  HXT(J)«IEXT(J+1) 

IEKP(IC)-KN 
GO  TO  100 


AD-A178  SSS 
UNCLASSIFIED 
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370  IF(JGBK3E.GT.O)  GO  ID  100 
C 

C  CALCULATION  OF  THE  OD&FICianS  OF  THE  BEST 
APRCKIMAT1GN 

C  USING  THE  INVERSE  DISCRETE  FOURIER  TRANSFORM 
C 

400  OONTINUE 
NM1-NFCNS-1 
FSEK1.0B-6 
GTQHM3RID(1) 

X(IBZ>— 2.0 

CN«2*NFCNS-1 

DHjm.O/QI 

LfI 

KKK4 

IF  (EDGE  (1)  .EQ.0.0.AND.EXX3E(2*EBAND6)  .R3.5)  KKK-1 
IF  (NEONS.  LE.  3)  KKK-1 
IFCKKK.  H3.1)  GO  TO  405 
DrEMP-DCDS(PI2*GRID(l) ) 

DNU!MXX)S(PI2*GRID(EGRID) ) 

AA-2 . 0/  (DTEMP-ENUM) 

BB—  (ETEMPfENUM)  /  (DTEMP-ENUM) 

405  CONTINUE 

DO  430  J-l,  NEONS 
ETXJ-1)  *DELF 
XMXDS(PI2*ET) 

IFOCKK.R3.1)  GO  TD  410 
XTV(XT-BB)/AA 
PRADOS  (XT)  /PI2 
410  XEWC(L) 

IF(XT.GT.XE)  GO  TD  420 
IF(  (XEMCT)  .I/T.FSH)  GO  TD  415 
I>Lfl 
GO  TD  410 
415  A(J)-YCL) 

GO  TD  425 

420  IF ( (XT-XE)  .IT. FSH)  GO  TD  415 
GRIDCD^T 
A(J)-GEE<1,IB) 

425  OONTINUE 

IF(L.GT.l)  LfL-1 
430  OONTINUE 

GRID(l)  -<3TBMP 

D0Q*»PI2/CN 

DO  510  >1,  NEONS 

DTEMPbO.O 

CNU^(J-1)*DDEN 

IF(JKL.LT.l)  GO  TD  505 

DO  500  K-i,IKL 

500  DTEKP»DTEMPfA(R+l)  *DOOS(ENUM*K) 

505  DTE»k>2.0«CTEMPfA(l) 

510  ALINA  (J) -CTEMP 
DO  550  J-2,  NEONS 


550  ALPHA(J)  -2*ALPHA(  J)  /CM 
ALPHA(l)  -ALPHA  (1)  /C N 


IF(KKK.Ba.l)  GO  TO  545 
P(l)  -2 . 0  *ALPHA  (NPCNS)  *BB+-ALPHA(NM1) 

P(2)  -2 . 0  *AA*ALPHA  ( NPCNS) 

Q(l)  -ALPHA  (NPCNS-2)  -ALPHA  (NPCNS) 

DO  540  J-2,l«L 
IP(J.UT.MCL)  GO  TO  515 
AAp0.5*AA 
BB-0.5*BB 
515  CONTINUE 
P(J+1)^).0 
DO  520  K-l/J 
AGO  «P(K) 

520  P(K)«2.0*BB*A(K) 

P(2)-P(2)4A(1)  *2.0*AA 
JMl-J-1 

DO  525  K-1/JM1 
525  P(K)  -P(K)  -K)(K)  +AA*A(K+1) 

JP1-J+1 

DO  530  K«3,JP1 
530  P(K)-P(K)4AA*A(K-1) 

IF(J.BQ.mi)  GO  TO  540 
DO  535  K-l/J 
535  QUO— A(K) 

Q(l)  -Q(l)  +ALPHA(NFCNS-1-J) 

540  CONTINUE 

DO  543  J-l, NPCNS 
543  ALPHA( J)  -P( J) 

545  (ENTINUE 

IF  (NPCNS.  GT.  3)  REHJFN 
ALPHA(NPCNSfl)  *0.0 
ALPHA  (NFCN9+-2)  *0 .0 
RETUBN 
END 
C 
C 
C 

DOUBLE  PRECISION  FUNCTION  D(K,N,M) 

C 

C  FUNCTION  TO  CALCULATE  OHE  LAGRANGE  INTERPOLATION 
C  ODEFFTCXQnS  FOR  USE  IN  THE  FUNCTION  GEE. 

C 

COMON 

PI2fAD,  DEV,  X,Y,  GRID,  EES,  WT,  ALPHA,  IEXT,  NPCNS,  N3RID 

DIMENSION  IEXT(514)  ,AD<514)  ,ALPHA(514)  ,X(514> ,  Y(514) 
DIM9ISIGN  EES (8224)  ,  GRID  (8224)  ,WT(8224) 

DOUBLE  PRECISION  AD,EEV,X,Y 
DOUBLE  PRECISION  Q 
DOUBLE  PRECISION  PT2 
D-1.0 
Q-X(K) 

DO  3  If1,M 
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DO  2  J-L,N,M 
IF (J-K) 1,2,1 

1  D-2.0«D*(Q-X<J)) 

2  CONTINUE 

3  OONTINUE 
D-1.0/D 
RETURN 
END 

C 

C 

c 

c 

DOUBLE  PRECISION  FUNCTION  GEE(K,N) 

C 

C  FUNCTION  TO  EYALQATO  THE  FRffiUBJCY  RESPONSE  USING  THE 
C  LAGRANGE  INTERPOLATION  FORHJLA  IN  THE  BARYCTNTRIC  FORM 
C 

COMMON 

H2 ,  AD,  DEV,  X,  Y,  GRID,  DES ,  WT,  ALKA,  TEXT,  NFCNS,  NGRID 

DDBJSICN  IEXr(514)  ,AD(514)  ,ALPHA(514)  ,X(514)  ,Y(514) 
DIMNSIGN  EES (8224)  ,GRID(8224)  ,WT(8224) 

DOUBLE  PRECISION  P,  C,D,XF 
DOUBLE  PRECISION  P32 
DOUBLE  PRECISION  AD,EEV,X,Y 
IMJ.O 
XF-GRID(K) 

XF-DO0S(PI2*XF> 

IX).  0 

DO  1  J»1,N 
OXF-X(J) 

OAD(J)/C 

DHDfC 

1  P»PfC*Y(J) 

GEE^P/D 

RETURN 

END 

C 

C 

C 

C 

SUBROUTINE  OUCH 
PRINT  1 

1  FORMAT  ('  ***********  FAILURE  TO  OOWERGE 

************  i  / 

1  'OERCBABLE  CAUSE  IS  MACHINE  RQUICING  ERROR1/ 

2  'OTHE  IMPULSE  RESPONSE  MAY  BE  OORRECT'/ 

3  _  'OGBBCK  WITH  A  FREQUENCY  RESPONSE' ) 

RETURN 

BID 


•WaVAvKsV,'' ‘.s' 


APPENDIX  B 

ZERO  EXTRACTION  PROGRAM 
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RPOLY  ROUTINE  BY  M.A.  JENKINS 
ACM  TOMS  VCL  1  NO  2  JUNE  1975 

SUBROUTINE  RPOLY 

FINDS  IHE  ZEROS  OP  A  REAL  POLYNOMIAL 
OP  -  DOUBLE  ERECTS  ION  VECTOR  OF  CDETFICTENTS  IN 
ORDER  OF  DECREASING  POWERS. 

DEGREE  -  INTEGER  DEGREE  OF  IHE  POLYNOMIAL. 

ZERCR,  ZEROI  -  OUTHJT  DOUBLE  IRBOSICN  VECTORS  OF 
REAL  AND  IMAGINARY  PARTS  OF  IHE 
ZEROS. 

FAIL  -  OUTHJT  LOGICAL  PARAMETER/  HUE  ONLY  IF 
LEADING  COEFFICIENT  IS  ZERO  OR  IF  RPOLY 
HAS  FOOT®  FEWER  IHAN  DEGREE  ZEROS. 

IN  IHE  LATTER  CASE/  DEGREE  IS  RESET  TO 
IHE  NUfBER  OF  ZEROS  POUND. 

TO  CHANGE  IHE  SIZE  OF  POLYNOMIALS  WHICH  CAN  BE 
SOLVED,  RESET  IHE  DIMENSIONS  OF  IHE  ARRAYS  IN  IHE 
COMMON  AREA  AND  IN  IHE  FOLLOWING  DECLARATIONS 
FOR  SCALING/  BOOTES  AND  ERROR  CALCULATIONS.  ALL 
CALCULATIONS  FOR  IHE  ITERATIONS  ARE  DONE  IN 
DOUBLE  PRECISION. 

ODWCN/FILE/  AMP (40 96)  /PHASE (40 96)  ,NFET,  ITYPE 
QOMMON/DAT/  FO,TFLO,T5HI,NUM 
OOMCN  /RPOLY/  P,  QP,  K,  QK/  SVK,  SR/  SI,  D, 

1  V,  A,  B,  C,  D,  AL,  A2,  A3,  A£,  A7,  E,  F,  G, 

2  H,  SZR,  SZI,  12 R/  121,  ETA,  ARE,  MRE,  N,  IN 
DOUBLE  PRECISION  P(1024) ,  QP(1024) ,  KC1024) , 

1  OR  (1024) ,  arc  (1024),  SR,  SI,  U,  V,  A,  B,  C,  D, 

2  Al,  A2,  A3,  A6,  A7,  E,  F,  G,  H,  SZR/  SZI, 

3  12  R,  12 1 
REAL  ETA,  ARE/  IRE 
INTEGER  N,  NN 

QP  IS  EfflBlSIONED  TO  4096  SINCE  IT  TAKES  ITS  ARRAY 
FROM  SAWCAD'S  AMP  ARRAY. 

DOUBLE  IRBCISION  CE>(4096) ,  TEMP(1024) , 

1  ZERCRQ024)/  ZEROK1024) ,  T,  AA,  BB,  OC,  DABS, 

2  FACTOR 

REAL  FT (1024) ,  ID,  MAX,  MIN,  XX,  YY,  CDSR, 

1  SINR,  XXX,  X/  SC/  BND,  XM,  FF,  DF,  DX,  INFIN, 

2  SMALNO,  BASE 

INTEGER  DEGREE/  COT,  IE,  I,  J,  JJ,  IKL 
LOGICAL  FAIL,  ZERCK 

THE  FOLLOWING  STATEMEOTS  SET  MACHINE  CONSTANTS  USED 
IN  IHE  VARIOUS  PARTS  OF  IHE  PROGRAM.  IHE  IGANING  OF  IHE 
FOUR  OONSTAtnS  ARE... 

ETA  IHE  IRXIMJM  RELATIVE  REPRESENTATION  ERROR 
WHICH  CAN  BE  DESCRIBED  AS  IHE  SMALLEST 
POSITIVE  FLOATING  POINT  NUM3ER  SUCH  THAT 


onojo  non  non  oooooo 
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NUM3ER  IF  THE  EXPONENT  RANGE  DIFFERS  IN  SINGLE 
AND  DOUBLE  ERECISION  THEN  SMALNO  AND  INFIN 
SHOULD  INDICATE  THE  SMALLER  RANGE. 

BASE  THE  BASE  OF  THE  FLOATING  POINT  NUMBER 
SYSTEM  USED. 

THE  VALUES  BELCH  CORRESPOND  TO  THE  VAX  11-750 
BASB*2. 

ETEBf1.387779E}-17 
INFINb1.7E38 
SMALN>5.9EH39 

ARE  AND  ERE  REFER  TO  THE  UNIT  ERROR  IN  +  AND  * 
RESPECTIVELY.  THEY  ARE  ASSUMED  TO  BE  THE  SAME  AS 
ETA 

ARE^ETA 
MlErElA 
LOSMALNQ/ETA 

INITIALIZE  ARRAYS 

IF  (ITYPE.EQ.l)  THEN 

HUNT  *,  'THIS  IS  A  FREQUENCY  FILE  I ' 

HUNT  *,  'EXPECTED  AN  IMPULSE  RESPONSE  FILE' 
RETURN 

END  IF 
D0GREE>NUM-1 
DO  233  1=1,  DEGREE*- 1 
QP(I)-AMP(I) 

AMP(I)4).0D0 
IHASE(I)«0.0D0 
ZERGR(I)*O.ODO 
ZEROKDm.ODO 
33  CONTINUE 

INITIALIZATION  OF  CONSTANTS  FOR  SHIFT  ROTATION 

XX-0 .70710678 
YY--XX 

OQSI*-.069756474 
SUB*. 997 56 405 
FAIL*.  FALSE. 

{^DEGREE 
NN-Nfl 

C  ALGORITHM  FAILS  IF  THE  LEADING  COEFFICIENT  IS  ZERO. 

IF  (OP(l)  .NE.O.DO)  GO  TO  10 
FAIL-.  TRUE. 

EBGRED«0 

HUNT  *,  'LEADING  CDEFFICIINT  IS  ZERO  —  NOT  ALLOWED' 
RETURN 

C  REMOVE  THE  ZEROS  AT  THE  ORIGIN  IF  ANY 
10  IF  (OP  (NO  .NE.0.0D0)  GO  TO  20 
J-DEGREEhNtl 
ZERCR(J)«O.ODO 
anp(j)^).0d0 
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ZHOI(J)  -O.ODO 
phase  (j)m.OdO 
tt*=MKL 
B-B-l 
GO  ID  10 

C  MAKE  A  COPY  OF  THE  ODEFFICXENIS 

20  DO  30  I-1,NN 
P(I)*X>(I) 

30  OONTINDE 

C  START  THE  ALGORITHM  TOR  ONE  ZERO 

40  IF(N.CT.2)GO  TO  60 
IF  (N.LT.1)  RETURN 

C  CALCULATE  THE  FINAL  ZERO  CR  FAIR  OF  ZEROS. 
IF(N.BQ.2)00  TO  50 
ZBVXKDBGREE)  — P(2)  /P(l) 

ZH8QI  (DEGREE}  -O.ODO 
AMP  (DEGREE)  «ZER0R  (DEGREE) 

PHASE  (DEGREE)  «ZERDI  (DEGREE) 

RETURN 

50  CALL  QQAD(P(1) ,  P(2) ,  P(3)  ,  ZERCR(DBGREE-l) , 

1  ZEROI(DBGREB-l) ,  ZERCR  (DEGREE) ,  ZEROI  (DEGREE) ) 
amp(degree-l)  -zeror  (degree-1) 
phase  (degree-1)  -zerol  (degree-1) 

AMP  (DEGREE)  «ZERCR(DBGREE) 

PHASE  (DEGREE)  “ZESOI  (nBEREEl 
RETURN 

C  FIND  THE  LARGEST  AND  SMALLEST  MODULI  OF  CDEFFTCIQnS 

60  MAX-O. 


MIN-INETN 
DO  70  I-1,NN 
X-ABS(SNGL(P(I))) 

IFOC.GT.MAX)  MAX-X 
IF(X.NE.O.  .AND.  X.  IT.  MIN)  MIN-X 
70  CONTINUE 

C  SCALE  IF  THERE  ARE  LARGE  OR  V&U  SMALL  (DHTICXENTS 
C  CDMRJTES  A  SCALE  FACTOR  TO  MULTIFLY  TEE 
C  COEFFICIENTS  OF  THE  POLYNOMIAL.  THE  SCALING  IS  DONE 
C  TO  AVOID  OVERFLOW  AND  TO  AVOID  UNCETBCTED  UNDERFLOW 


110  DO  120  >l,m 

PT(I)  W©SOCL(P(I) ) ) 

120  CONTINUE 

ETOW)  — ET(NN) 

C  COMPUTE  UPPER  ESTIMATE  OF  BOUM) 

X-EXP  ( (ALOG  ( -FT  (NN) )  -ALOG  (PT(1) ))  /FLOAT  (N) ) 
IF(PT(N)  .Ha.OJQO  TO  130 

C  IF  NOTION  STEP  AT  THE  ORIGIN  IS  BETTER/  USE  IT! 
X^-PT(M0/PT(N) 

IF  (XM.  LT.  X)  X-XM 

C  CHOP  THE  INTERVAL  (0,X)  UNTIL  FF  .LE.  0 
130  X»MC*.l 
FF-PT(l) 

DO  140  1-2, NN 
FF-FE*XHHT(I) 

140  CONTINUE 

IFtFF.LE.O .) 00  TO  150 

y-YH 

GO  TO  130 
150  DX-X 

C  DO  NWTEN  ITERATION  UNTIL  X  CONVERSES  TO  TWO 
C  DECIMAL  PLACES 

160  IF(ABS(DX/X)  .LE.  .005)00  TO  180 
FF-ET(l) 

DF-FF 

DO  170  1-2, N 
FF*FF*X+FT<I) 

DF"OF*X+FF 
170  CONTINUE 

FF-FF*X+PT(NN) 

QX-FF/DF 
X»*-OX 
GO  TO  160 
180  BND-X 

C  COMPUTE  THE  DERIVATIVE  AS  THE  INITIAL  K  POLYNOMIAL 
C  AND  DO  5  STEPS  WITH  NO  SHIFT 
NM1-N-1 
DO  190  1-2 ,N 

K(I)«FLQAT(m-I)*P(I)/FLQAT(N) 

190  OONTINUE 
K(l)-P(l) 

AAfP(NN) 

B&*P(N) 

ZERCK«F(N)  .H3.0.D0 
DO  230  J>lf5 
OOK(N) 

IF(ZEROK)GO  TO  210 


C  USE  SCALED  FORM  OF  RECURRENCE  IF  VALUE  OF  K  AT  0  IS 
C  NONZERO 
TV-AA/CC 
DO  200  1-1, NM1 
J-NN-I 

K(J)-T*K(J-1)+P(J) 


■ 
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200  CONTINUE 
K(l)-P(l) 

ZER0K4ABS(K (N) )  .LE.DftBS(BB)  *ETA*10. 

GO  TO  230 

C  USE  THE  UNSCALED  FORM  OF  RECURRENCE 
210  DO  220  I«1,MCL 
J*NN-I 
K(J)  “K(J-l) 

220  CONTINUE 
K  (1)0.  DO 
ZERQfWK(N)  .EQ.O.DO 
230  CONTINUE 

C  SAVE  K  FOR  RESTARTS  WITH  N0I  SHIFTS 
DO  240  1-1 ,N 
TEMP(I)»K(I) 

240  CONTINUE 

C  LOOP  TO  SELECT  THE  QUADRATIC  CORRESPONDING  TO  EACH 
C  NEW  SHIFT 

DO  280  anvi,20 

C  QUADRATIC  RESPONDS  TO  A  DOUBLE  SHIFT  TO  A 
C  NON- REAL  POINT  AND  ITS  COMPLEX  CONJUGATE.  THE  POINT 
C  HAS  MODULUS  BND  AND  AMR. DUDE  ROTATED  BY  94  DEGREES 
C  FROM  THE  IREVIDUS  SHIFT 
XXXKDSR*XX-SINR*YY 
YY«SINR*XX-KDSR*YY 

xx-xxx 

SI*BND*XX 

SI«BND*YY 

D—2.0D0*SR 

V-BND 

C  SECOND  STAGE  CALCULATION,  FIXED  QUADRATIC 
CALL  FXSHFR(20*CNT,FE) 

IF(NZ.HQ.O)GO  TO  260 

C  THE  SECOND  STAGE  JUMPS  DIRECTLY  TO  ONE  OF  THE  THIRD 
C  STAGE  ITERATIONS  MID  RETURNS  HERE  IF  SUCCESSFUL. 

C  DEFLATE  THE  POLYNOMIAL,  STORE  THE  ZERO  CR  ZEROS  AND 
C  RETURN  TO  THE  MAIN  ALGORITHM. 

J-DBGREEKNfl 
ZERCR(J)  -SZR 
AMP(J)  -SZR 
ZEROI(J)  -SZI 
EHASE(J)  -SZI 
NfrNN-lE 
N»NN-1 

DO  250  I*1,NN 
P(I)KJP<I> 

250  CONTINUE 

IF(FE.EQ.l)  GO  TO  40 
ZERCR(Jfl) -I2R 
ZEROI  (Jfl)  -12 1 
AMP(J+1)  “I2R 
PHASE  (Jfl) -12 1 
GO  TO  40 
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C  IF  THE  ITERATION  IS  UNSUCCESSFUL  MOTHER  QUAERATIC 
C  IS  CHOSEN  AFTER  RESTORING  K 
260  DO  270  I«1,N 
K(I)-TSMP(I) 

270  CONTINUE 
280  CONTINUE 

C  RETURN  WITH  FAILURE  IF  NO  CONVERGENCE  WITH  20 
C  SHOTS. 

FAIL- .  TRUE. 

DBGREEWDEJGREEhN 

RETURN 

END 

SUBROUTINE  FXSHER(L2,  EE) 

C  OOMRJTES  UP  TO  12  FIXED  SHIFT  K-POLYNCMIALS, 

C  TESTING  FOR  CONVERGENCE  IN  THE  LINEAR  CR  QUAERATIC 
C  CASE.  INITIATES  CNE  OF  THE  VARIABLE  SHIFT 
C  ITERATIONS  AND  RETURNS  WITH  THE  NUMBER  OF  ZEROS 
C  FOUND. 

C  L2  -  LIMIT  OF  FIXED  SHIFT  STEPS. 

C  NZ  -  NUMBER  OF  ZEROS  FOUND. 

COMO*  /RFCLY/  P,  QP,  K,  QR,  SVK,  SR,  SI,  U, 

1  V,  A,  B,  C,  D,  Al,  A2,  A3,  A£,  A7,  E,  F,  G, 

2  H,  SZR,  SZI,  12  R,  121,  ETA,  ARE,  MLE,  N,  NN 
DOUBLE  PRECISION  PQ024) ,  QPQ024) ,  KQ024) , 

1  QR (1024) ,  SVK(1024) ,  SR,  SI,  U,  V,  A,  B,  C,  D, 

2  AL,  A2,  A3,  A6,  A7,  E,  F,  G,  H,  SZR,  SZI, 

3  12  R,  12 1 

REAL  ETA,  ARE,  MCE 
INTEGER  N,  NN 

DOUBLE  ERBCISICN  SVU,  SW,  UI,  VI,  S 
REAL  BETAS,  BET3BV,  OSS,  CW,  SS,  W,  TS,  TV, 

1  GTS,  OTV,  TW,  TSS 

INTEGER  12,  IB,  TYPE,  I,  J,  IFLAG 

LOGICAL  VPASS,  SPASS,  VTRY,  STRY 

IB-0 

BEEBV-.25 
BETA9-.25 
069“  SR 

aw-v 

C  EVALUATE  POLYNOMIAL  BY  SYNTHETIC  DIVISION 
GALL  QUADSD(NN,U,V,P,QP,A,B) 

CALL  CALGSC(TYPE) 

DO  80  J-1,12 

C  CALCULATE  NEXT  K  POLYNOMIAL  AND  ESTIMATE  V 
CALL  NEXTR(TYPE) 

CALL  CALCSC(TYPE) 

CALL  N9VEST(TYPE,  UI,VI) 

W*VI 

C  ESTIMATE  S 
S9>0. 

IF(K(N)  .NE.O.DO)  S>-P(N©/K(N> 

TV-1. 
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IF(J.B2.1  .OR.  TXPE.BQ.3)  GO  TO  70 
C  OOMTOTE  RELATIVE  MEASURES  OF  CONVERGENCE  OF  S  AND  V 
C  SEQUENCES 

IF(W.  NE.O.)  TV*ABS(  (W-OW)  /W) 
IF(SS.NE.O.)TS*ABS( (SS-OSS)/SS) 

C  IF  INCREASING,  MJLTIH.Y  TWO  MOST  RECENT 
C  CONVERGENCE  MEASURES 
TW*1. 

IF  (TV.  LT.  OTV)  TW=TV*OTV 
TSS-1. 

IF(TS.LT.CT5)TS&'TS*0'IS 
C  COMPARE  WITH  CONVERGENCE  CRITERIA 
Tr2AS&«TW.  LT.  BETSV 
SPAS9*T5S.  LT.  BETAS 
IF  (.NOT.  (SPASS  .OR.  VEASS) )  GO  TO  70 
C  AT  LEAST  ONE  SEQUENCE  HAS  PASSED  THE  CONVERGENCE 
C  TEST.  STORE  VARIABLES  BEFORE  ITERATING. 

SVl^U 
SW*V 

DO  10  1*1, N 
SVK(I)-«(I) 

10  CONTINUE 
9-SS 

C  CHOOSE  ITERATION  ACCORDING  TO  THE  FASTEST 


C  CONVERGING  SEQUENCE. 

VTRY-.  FALSE. 

STRY* .  FALSE. 

IF  (SPASS  .AND.  ( (.NOT.  VEASS)  .CR. 

1  TSS.LT.TW) )  GO  TO  40 
20  CALL  QUADtT(UI,VI,IC) 

IF(EC.GT.O)  RETURN 

C  QUADRATIC  ITERATION  HAS  FAILED.  FLAG  THAT  IT  HAS 
C  BEEN  TRIED  AND  DECREASE  THE  CONVERGENCE  CRITERION. 
VTRYfc.TOTE. 

BETHV*BETAV*.25 

C  TRY  LINEAR  ITERATION  IF  IT  HAS  NOT  BEEN  TRIED  AND 
C  THE  S  SEQUENCE  IS  CONVERGING. 

IF(STRY  .OR.  ( .NOT. SPASS) )  GO  TO  50 

DO  30  1*1, N 

R(I)-SVK(I) 

30  CONTINUE 

40  CALL  REAL  IT  (S,  EC,  IFLAG) 

IF(NZ.GT.O)  RETUFN 

C  LINEAR  ITERATION  HAS  FAILED.  FUG  THAT  IT  HAS  BEEN 
C  TRIED  AND  DECREASE  THE  CONVERGENCE  CRITERION. 

STRY“  .TRUE. 

BETAS- BETAS *.25 
IF  (IFLAG.  H3. 0)  GO  TO  50 

C  IF  LINEAR  ITERATION  SIGNALS  AN  ALMOST  DOUBLE  REAL 
C  ZERO  ATTEMPT  QUADRATIC  ITERATION. 

UI— (SfS) 

VI-S*S 
GO  TO  20 
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C  RESTORE  VARIABLES 
50  OSVO 

V=SW 

DO  60  1*1, N 
K(I)*S7K(I) 

60  CONTINUE 

C  my  QUADRATIC  ITERATION  IF  IT  HAS  NOT  BEEN  miED 
C  AND  THE  V  SEQUENCE  IS  CONVERGING. 

IF(VPASS  .AND.  (.NOT.VTOY))  GO  TO  20 
C  REXDMHHE  QP  AND  SCALAR  VALUES  TO  CONTINUE  THE 
C  SECOND  STAGE. 

CALL  QUADSD(NN,  U,  V,  P,  QP,  A,  B) 

CALL  CALCSC(TYPE) 

70  OW=W 

OSS-SS 
OTV*TV 
OIS*TS 

80  CONTINUE 

RETOIW 
END 

SUBRCUTINE  QUADIT(OU,W,JE) 

C  VARIABLE- SHIFT  K-POLYNCMIAL  ITERATICN  FOR  A 
C  QUADRATIC  FACTOR  CONVERGES  CNLY  IF  THE  ZEROS  ARE 
C  EQUIMDEULAR  OR  NEARLY  SO. 

C  UU,  W  -  COEFFICIENTS  OF  STORTING  QUADRATIC 
C  NZ  -  NUEBER  OF  ZEROS  FOUND 

COMMON  /RPOLY/  P,  QP,  K,  QK,  SVK,  SR,  SI,  D, 

1  V,  A,  B,  C,  D,  AL,  A2,  A3,  A6,  A7,  E,  F,  G, 

2  H,  SZR,  SZI,  IZR,  IZI,  ETA,  ARE,  «E,  N,  NN 
DOUBLE  HiEOSION  P(1024) ,  QP(1024) ,  K(1024) , 

1  QR (1024) ,  SVK(1024),  SR,  SI,  D,  V,  A,  B,  C,  D, 

2  Al,  A2,  A3,  A 6,  A7,  E,  F,  G,  B,  SZR,  SZI, 

3  IZR,  IZI 

REAL  ETA,  ARE,  HIE 
INTEGER  N,  NN 

DOUBLE  IREXUSICN  DI,  VI,  DU,  W,  DABS 
REAL  (6,  MP,  CMP,  EE,  RELSTP,  T,  ZM 
INTEGER  m,  TYPE,  I,  J 
LOGICAL  TRIED 
NZ*0 

TRIE3>>.  FALSE. 

O-UU 

V-W 

J*0 

C  MAIN  LOOP 

10  CALL  QUADd.DO,  U,  V,  SZR,  SZI,  IZR,  IZI) 

C  RETUFN  IF  ROOTS  OF  IHE  QUADRATIC  ARE  REAL  AND  NOT 
C  CLOSE  TO  MULTIPLE  OR  NEARLY  EQUAL  AND  OF  OPPOSITE 
C  SIGN. 

IF  (DABS  (DABS (SZR)  -TABS(IZR) )  .GT.  .01D0* 

1  DABS(IZR) )  RETUR] 

C  EVALUATE  POLYNOMIAL  BY  QUADRATIC  SYNTHETIC  DIVISION, 
CALL  QUADSD(NN,  U,  V,  P,  QP,  A,  B) 
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MP»DABS(A-SZR*B)  «ABS(SZI«B) 

C  CDMHJTE  A  RIGCPOJS  BGU1D  ON  THE  RCU1DIN3  H®CR  IN 
C  EVALUATING  P 

Zft*SQFT(ABS(  £NGL(  V) ) ) 

EE>2.*ABS(SNGL(QP(1) ) ) 

tv-szr*b 
CO  20  1-2, N 

E2>EE*ZMfABS(£NGL(QP(I) ) ) 

20  CONTINUE 

EEXEE*ZMfABS(SN3L(A)  +T) 

ED-(5.*MREf4.*ARE)  *EE“(5.«MREf2.*ARE)  * 

1  (ABS(£NGL(A)  +D  +ABS(SNGL(B) )  *ZM)  + 

2  2.*ARE*ABS(T) 

C  ITERATION  HAS  CONVERSED  SUFFIGENILY  IF  THE 
C  POLYNOMIAL  VALUE  IS  LESS  THAN  20  TDES  IE IS  BOJND 
IF(MP.GT.20.*EE)GO  TO  30 
EC-2 
RETURN 

30  J-Jfl  _ 

C  STOP  ITERATION  AFTER  20  STEPS 
IF  (J.GT.  20)  RETURN 
IF (J.  IT. 2)  GO  TO  50 

IF(RaSTP.GT..01  .OR.  MP.LT.CMP  .OR.  TRIED) 
1  GO  TO  50 

C  ACLUSTER  APPEARS  TO  BE  STALLING  THE  CONVERGENCE. 

C  FIVE  FIXED  SHIFT  STEPS  ARE  TAKEN  WITH  A  U,V  CLOSE 
C  TO  THE  CLUSTER. 

IF  (RELSIP.LT.  ETA)  RELSTWETA 
RELCTP-  32  RT  (RELSTP) 

(>»U-U*RaSTP 

V-V4-V*RELSTP 

CALL  QUAD6D(NN,U,V,P,QP,A,B) 

DO  40  1-1,5 
CALL  CALCSC(TYPE) 

CALL  NEKCK(TYPE) 

40  OQNTINUE 

TRIEEN.THJE. 

J-0 

50  CMP-MP 

C  CALCULATE  NEXT  R  POLYNOMIAL  AND  NW  U  AND  V 
CALL  CALCSC(TYFE) 

CALL  NEXTK(TYPE) 

CALL  CALCSC(TYIE) 

CALL  NWEST(TYPE,  UI,  VI) 

C  IF  VI  IS  ZERO  THE  ITERATION  IS  NOT  CONVERGING. 

IF (VI. HQ. 0. DO)  RETURN 
RELSTM»BS(  (VI-V)  /VI) 

U»UI 
V-VI 
GO  TO  10 
END 

SUBROUTINE  REALIT(SSS,  W,  HUG) 

C  VARIABLE  SHIFT  H  POLYNOMIAL  ITERATION  FOR  A  REAL 
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C  ZERO.  _ 

C  SSS  -  STARTING  ITERATE 
C  NZ  -  NUEBER  OP  ZERO  FOUID 

C  IELAG-  FLAG  TO  INDICATE  A  PAIR  OF  ZEROS  NEAR  REAL 
C  AXIS. 

OOMM2N  /RFOLY/  P,  QP,  K,  QR,  SVK,  SR,  SI,  U, 

1  V,  A,  B,  C,  D,  Al,  A2,  A3,  A£,  A7,  E,  P,  G, 

2  H,  SZR,  SZI,  12  R,  IZI,  ETA,  ARE,  ERE,  N,  W 
DOCBLE  PRECISION  PC1024) ,  QPQ024)  ,  KC1024) , 

1  QR (1024),  9/KQ024) ,  SR,  SI,  0,  V,  A,  B,  C,  D, 

2  Al,  A2,  A3,  A6,  A7,  E,  P,  G,  H,  SZR,  SZI, 

3  12  R,  IZI 

REAL  ETA,  ARE,  ERE 
INTEGER  N,  NN 

DOUBLE  PRECISION  IV,  XV,  T,  S,  SSS,  DABS 
REAL  (6,  MP,  OKP,  EE 
INTEGER  NZ,  IFLAG,  I,  J,  NKL 
MCL-N-1 
M«0 
&.SSS 
IFLAOO 

JmO 

C  MAIN  LOOP 
10  PV-P(l) 

C  EVALUATE  P  AT  S 
QP(1)-FV 
DO  20  1-2,  Mf 
FV-IV*S+P(I) 

QP(I) -IV 
20  OONTINQE 

MP^DABS(PV) 

C  COMPUTE  A  RIGOROUS  BOOED  ON  THE  ERROR  IN  EVALUATING 
C  P 

MS-DABS  (S) 

EEXMRE/ (AREMRE) )  *ABS(£NGL(QP(1)  ) ) 

DO  30  1-2, MI 

EE^EE*MS+ABS(£N3L(QP(I) ) ) 

30  OONTINUE 

C  ITERATION  HAS  OONVERGH)  SUFFICIENTLY  IF  THE 
C  POLYNOMIAL  VALUE  25  LESS  THAN  20  TIMES  THIS  BODED. 
IF(MP.GT.20.*(  (AREMRE)*EB-MRE*MP))QO  TD  40 
NZ-1 
SZRiS 
SZI-O.DO 
RETURN 

40  JOfl  _ 

C  STOP  ITERATION  AFTER  10  STEPS 
IF(J.GT.IO)  REHJMI 
IF(J.LT.2)  GO  TO  50 

IF  (DABS  (T)  .GT.  .001*DABS(S-T)  .OR.  MP.LE.OMP) 

1  GO  TO  50 

C  A  CLUSTER  OP  ZEROS  NEAR  THE  REM,  AXIS  HAS  BEEN 
C  ENCOUNTERED.  RETURN  WITH  IFLAG  SET  TD  INITIATE 
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C  QUADRATIC  ITERATION. 

DUG-1 

SS£>S 

RETURN 

C  RETURN  IF  THE  POLYNOMIAL  VALUE  BAS  INCREASED 
C  SIGNIFICANTLY. 

50  CME*MP 

C  COMPUTE  T,  THE  NEXT  POLYNOMIAL,  AND  THE  NBV  HERAT 
KV-K(l) 

QK(1) *KV 
DO  60  1-2 /N 
FV«KV*SfK(I) 

QK(I)-RV 
60  OONTINUE 

IF  (DABS  (KV)  .LE.DABS(K(N))*10.*ETA)GO  TO  80 
C  USE  THE  SCALED  FORM  OF  THE  RECURRENCE  IF  THE  VALUE 
C  OF  K  AT  S  IS  NONZERO 
TV-IV/KV 
K(1)KJP(1) 

DO  70  1-2, N 
Kd)*r*QKd-l)4QP(I) 

70  OONTINUE 

GO  TO  100 
C  USE  SCALED  FORM 
80  K(1)-4.0D0 

DO  90  1-2, N 
K(I)^K(I-1) 

90  CONTINUE 

100  KV-K(l) 

DO  110  1-2, N 
KV-KV*SfK(I) 

UO  CONTINUE 

TVO.DO 

IF (DABS(KV)  .GT.  DABS (K  (N) )  *10 .  *ETA)  TV-PV/KV 

»SfT 

GO  TO  10 

END 

SUBROUTINE  CALCSC(TYPE) 

C  THIS  ROUTINE  CALCULATES  SCALAR  QUANTITIES  USED  TO 
C  COMPUTE  THE  NEXT  K  POLYNOMIAL  MD  NEW  ESTIMATES  OF 


C  THE  QUADRATIC  COEFFICIENTS. 

C  TYPE  -  INTEGER  VARIABLE  SET  HERE  INDICATING  HCN  THE 
C  CALCULATIONS  ARE  NORMALIZED  TO  fl/OID 

C  OVERFLOW. 

COMMON  /RFQLY/  P,  QP,  K,  QK,  SVK,  SR,  SI,  U, 

1  V,  A,  B,  C,  D,  Al,  A2,  A3,  A6,  A7,  E,  F,  G, 

2  H,  SZR,  SZI,  12  R,  12 1,  ETA,  ARE,  IRE,  N,  W 
DOUBLE  PRECISION  P(1024) ,  QP(1024)  ,  K(1024) , 

1  QM  (1024),  SVK(1024),  SR,  SI,  U,  V,  A,  B,  C,  D, 

2  Al,  A2,  A3,  A6,  A7,  E,  F,  G,  B,  SZR,  SZI, 

3  12  R,  121 

REAL  ETA,  ARE,  IRE 


n  n> 
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DOUBLE  PRECISION  DABS 
INTEGER  TYPE 

C  SYNTHETIC  DIVISION  OF  K  BY  THE  QUADRATIC  1,U,V 
CALL  QUADSD(N,  U,  V,  K,  QK,  C,  D) 

IF  (DABS  (C)  .GT.DABS(K(N) )  *100 .  *ETA)  QO  TO  10 
IF  (DABS  (D)  .GT.DABS(KCN-l) )  *100.*ETA)  GO  TO  10 
TYPE-3 

C  TYPE-3  INDICATES  THE  QUADRATIC  IS  ALMOST  A  FACTOR 
C  OF  K. 

RETURN 

10  IF  (DABS  (D)  .LT.DABS(C)  )GO  TO  20 

TYPE-2 

C  TYPE-2  INDICATES  THAT  ALL  FORMULAS  ARE  DIVIDED  BY  D 
E-A/D 
F-C/D 
G-U*B 
»«V*B 

A3-(AKJ)*DfH*(B/D) 

A1-B*F-A 
A7“(F-KD  *AfH 
RETURN 
0  TYP&'l 

TYPE-1  INDICATES  THAT  ALL  EORHJLAS  ARE  DIVIDED  BY  C. 
£>A/C 
F-O/C 
G-U*E 
B«V*B 

A3-A*Bf  (H/C+G)  *B 
A1-B-A*(D/C) 

A7-A+G*DfH*F 
RETURN 
END 

SUBROUTINE  NEXTK(TYPE) 

C  COMPUTES  THE  NEXT  K  POLYNOMIALS  USING  SCALARS 
C  COMPUTED  IN  CALCSC. 

COMMON  /RP0LY/  P,  QP,  K,  QK,  SVK,  SR,  SI,  U, 

1  V,  A,  B,  C,  D,  Al,  A2,  A3,  A6,  A7,  E,  F,  G, 

2  H,  SZR,  SZI,  LZR,  121,  ETA,  ARE,  HIE,  N,  HJ 
DOUBLE  PRECISION  P(1024) ,  QP(1024) ,  K(1024) , 

1  QR (1024) ,  SVK(1024) ,  SR,  SI,  D,  V,  A,  B,  C,  D, 

2  Al,  A2,  A3,  A6,  A7,  E,  F,  G,  B,  SZR,  SZI, 

3  LZR,  LZI 

REAL  ETA,  ARE,  IRE 
INTEGER  N,  NN 

DOUBLE  PRECISION  TEMP,  DABS 
INIHGER  TYPE 
IF  (TYPE.  53-3)  00  TO  40 
TSM^A 

IF  (TYPE. S3. 1)  TEMP-B 

IF  (DABS  (Al) .  GT.  DABS  (TEMP)  *ETA*10.)  GO  TO  20 
C  IF  Al  IS  NEARLY  ZERO  THEN  USE  A  SPECIAL  FORM  OF  THE 
C  RECURRENCE. 

K(1)4).D0 


104 


K(2)— A7*QP(1) 

DO  10  1-3, N 

K(I)WS*QK(I-2)-A7«QP(I-l) 

10  CONTINUE 

RETURN 

C  USE  SCALED  FORM  OF  THE  RECURRENCE. 

20  A7-A7/A1 

A3-A3/A1 
K(1)^JP(1) 

K(2)K)Pt2)-A7*QP(l) 

DO  30  1-3 ,N 

K(I)^A3«QK(I-2)  -A7«QP(I-1)  -tQP(I) 

30  CONTINUE 

RETURN 

C  USE  UNSCALED  FORM  OF  THE  RECURRENCE  IF  TYPE  IS  3 
40  K(1)«4).D0 

K(2)-O.DO 
DO  50  1-3 ,N 
K(I)»<*(I-2> 

50  CONTINUE 

RETURN 
END 

SUBROUTINE  NEWEST  (TYPE,  UU,  W) 

C  03MHJTES  NEW  ESTIMATES  OF  THE  QUADRATIC  COEFFICIENTS 
C  USHC  THE  SCALARS  CDMRJTED  IN  CALCSC. 

COMMON  /RPQLY/  P,  QP,  K,  CM,  £VK,  SR,  SI,  U, 

1  V,  A,  B,  C,  D,  AL,  A2,  A3,  A6,  A7,  E,  F,  G, 

2  H,  SZR,  SZI,  IZR,  12 1,  ETA,  ARE,  MRE,  N,  EN 
DOUBLE  EREdSIDN  P(1024) ,  QP(1024) ,  KC1024) , 

1  CM  (1024) ,  S/K(1024) ,  SR,  SI,  U,  V,  A,  B,  C,  D, 

2  Al,  A2,  A3,  A6,  A7,  E,  F,  G,  H,  SZR,  SZI, 

3  LZR,  IZI 

REAL  Eta,  ARE,  MRE 
INTBGHIN,  NN 

DOUBLE  ERBOSICN  A4 ,  A5  ,B1  ,R2  ,C1,C2,C3, 

1  04,  TEMP,  DU,  W 

INTEGER  TYPE  _ 

C  USE  FORMULAS  APHRCHUAIE  TO  SETTING  TYPE. 

IF  (TYPE. EQ. 3)  GO  TO  30 
IF  (TYPE.  EQ.  2)  GO  TO  10 
AA-AfrU*B+H*F 
A5«0(CH-V*F)  *D 
GO  TO  20 

10  A4-(ArK5)  *F+H 

A5-  (F+U)  *C+V*D 

C  EVALUATE  NEW  QUADRATIC  COEFFICIENTS. 

20  Bl— K(N)/P(N» 

B2— (K  (N-l)  +Bl*P  (N) )  /P  (NN) 

Cl«V*B2*Al 

C2»Bl*A7 

C3«B1*B1*A3 

C4-CL-C2-C3 

TTME*A5+B1*A4-C4 


IF  (TEMP.  03*0. DO)  GO  TO  30 
UtMJ-(U*(C3+C2)  +V*(Bl*Al+B2*A7)  )/TEMP 
W-V*  (1  .+C4/TEMP) 

RETURN 

C  IF  TYPE-3  TOE  QUADRATIC  IS  ZEROED. 

30  UOO.DO 
WO  .DO 
RETURN 
END 

SUBROUTINE  QUADSDCNN,  0,  V,  P,  Q,  A,  B) 

C  DIVIDES  P  BY  TOE  QUADRATIC  1,U,V  PLACING  TOE 
C  QUOTIENT  IN  Q  AND  TOE  REMAINDER  IN  A,B 

DOUBLE  PRECISION  P(NN) ,  Q(NN)  ,  0,  V,  A,  Bf  C 

INTEGER  I 

B-P(I) 

Q(1)«B 

AHP(2)-U«B 

Q(2)-A 

DO  10  I«3,EN 

OP(I)-U*A-V*B 

Q(I)-C 

B=A 

AfC 

10  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  QUAD  (A,  Bl,  C,  SR,  SI,  LR,  LI) 

C 

C  CALCULATE  TOE  ZENOS  OF  TOE  QUADRATIC  A*Z**2+B1*Z+C. 
C  TOE  QUADRATIC  FQRHJLA,  MODIFIED  TO  AVOID 
C  CVERFLCW,  IS  USB)  TO  FIND  TOE  LARGER  ZERO  IF  TOE 
C  ZEROS  ARE  REAL  AND  BOTH  ZEROS  ARE  COMPLEX. 

C  TOE  SMALLER  REAL  ZERO  IS  FOUND  DIRECTLY  FROM  TOE 
C  PRODUCT  OF  TOE  ZEROS  C/A. 

DOUBLE  PRECISION  A,  Bl,  C,  SR,  SI,  IR,  LI,  B, 
1  D,  E,  DABS,  DSQRT 

IFCA.NE.0.D0)  GO  TO  20 
SlH).D0 

IF(Bl.NE.O.DO)  Sft— C/Bl 
IA-O.DO 
10  SI-0. DO 

LI-0  .DO 
RETURN 

20  IF(C.NE.0.D0)  GO  TO  30 

SR-O.DO 
IA*-B1/A 
GO  TO  10 

C  COMPUTE  DISCRIMINATE  /VOIDING  OVERFLOW. 

30  M&/2.D0 

IF  (DABS  (B)  .LT.DABS(C) )  GO  TO  40 
E>»1  .DO-tA/B)  *(C/B) 

D-DSQRT(DABS(E) )  *DABS(B) 
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GO  TO  50 
40  E>=A 

IE'(C.LT.O.DO)  B»- A 
B«B*(E/DftBS(C))-E 
ENDSQRT(DABS(E) )  *DSQKT(DABS(C) ) 
50  IF(E.LT.O.DO)  GO  TO  60 

C  REAL  ZEROS. 

1F(B.  GE.O.DO)  D— D 
LR«{~BfD)/A 
Sl^O.DO 

IE'(LR.NE.O.DO)  SR-(C/LR)/A 
GO  TO  10 

C  COMPLEX  OCNJOGATE  ZEROS. 

60  Si^-a^A 

LR-SR 

SI«DABS(D/A) 

LI— SI 
RETDBN 
END 


APPENDIX  C 

OPTIMAL  COMBINATION  AND  RECONSTITUTION  PROGRAM 


i 

I 

I 

i 

) 


BE 


SUBPCUTINE  (DffiO 

c 

C  ISIS  ROUTINE  GENERATES  ALL  POSSIBLE,  NON  REPEATING 
C  OOWINATIDNS  OF  n  sameles  taken  k  at  a  time. 

C  ALGORITHM  ADAPTED  FROM 
C  APPLIED  COMBINATORIAL  MATHEMATICS 
C  RAGE  24 

C  POLYA  ET  AL  1964 

C  JOHN  WILEY  AND  SCNS,  INC,  NEW  YORK 

CDMCN/FILE/  AMP(4096) ,  PHASE  (40 96)  ,NFFT,  ITYPE 

(XXVQI/DAT/FO,TFLO,TFHI,NUM 

DOUBLE  PRECISION  ZER0R(1024) ,  ZEROK1024) 

DOUBLE  PRECISION  ffiEST0NE<1024) ,  IBESTTWO(1024) 
INTEGER  C(1024) ,  G(1024) ,  A(1024) 

INTEGER  I,K,N,T 

INTEGER  DEGREE,  TOTAL,  MONEBEST,  MWCBEST 


IF  (ITYEE.NE.0)  THEN 

HUNT  *,'***  NEED  TO  OBTAIN  ZEROS  FIRST 
RETURN 

END  IF 
DBGREE>NUM-1 
IF  (DEGREE.  GT.  35)  THEN 
HUNT  V  ' 

HUNT  *,'  ******  WARNING  ******* 

HUNT  *, 1  COMPUTATION  TIME  WILL  EXCEED  2 

HOURS1 

HUNT  *,'  * 

END  IF 

DO  289  1*1, DEGREE 
ZERGR(I)*AMP(I) 

ZEROKI)-PHASE(I) 

289  CONTINUE 
TOTAJ>0 
tW®GREE/2+2 

IF  ((ETGAT(N) -FLOAT (DEGREE) /2.0ED)  -NE.O)  THEN 
MM 

ZEROR(DBGREBfl)  -O.ODO 
ZEROI  (DBGREBfl)  -O.ODO 

END  IF 


C  INITIALIZE  ARRAY  C 
C 

DO  140  1*1, N 
C(I)*I 
140  CONTINUE 
C 

DO  5  K«0,N/2 

HUNT  *,K, '  AT  A  THE  FOR  HI 1 

HUNT  *,  TOTAL,  *  CD®  BUTTONS  TESTED  SO  EAR* 

IF  (K.E3.0)  GO  TO  1200 


109 


C  INITIALIZE  ID  BEGIN 
C 

mi 

A(l)  *1 

300  G(TX(A(T)) 

IP  (T.03.K)  GO  TO  1200 
A(TH)^A(T)+1 

IP  (A(Tfl)  .EQ.  (1+N) )  GO  TO  900 
mTfl 
GO  TO  300 
900  M-l 

IP  (T.BQ.O)  GO  TO  5 
GO  TO  1300 
1200  CALL 

POLYRBQQNC  ZERCR,  ZERDI,G,  DEGREE,  K,  ffiESTCNE, 

1  ffiESTIWO,  MDNEBEST,  MTWCBEST) 

TOTAL*  TOTAL+1 
IF  (K.TO.0)  GO  TO  5 
1300  A(T)-A(T)+1 

IP  (A(T)  .EQ.  (1+N) )  GO  TO  900 

GO  TO  300 
5  CONTINUE 

HUNT  *, '  ***  ALL  rTERATICNS  COMPLETE  ***' 

HUNT  *, 1  TOTAL  COMBINATIONS  ■  1 , TOTAL 

HUNT  *,  'BEST  DESIGN  FOUND: 1 
HUNT  ' 

HUNT  *, 'TRANSDUCER  1:' 

HUNT  • 

DO  888  Ol,M0NEBEST 

HUNT  *,  'H(*,L, ')  «MBESTONE(L) 

AMP(L)  =BBESTCNE(L) 

PHASE  (D  =>0.0 

888  CONTINUE 
itype*-i 

NUJMCNEBEST 
NFFMENEBEST 
KM).  5 

TET/>-(NCJM-l)/2.0 
TJHI*  (NUM-1)  /2.0 
CALL  WRITBO 
HUNT  ' 

HUNT  *,  'TRANSDUCER  2: ' 

HUNT  V  ' 

DO  889  I/»1,MIWCBEST 

HUNT  *,  'H(',L,')  -',ffiESTIWO(L) 

AMP(L)  -SBESTIWOdJ 
PHASE(L)  *0.0 

889  CONTINUE 
mro-i 
NUlMnWCBEST 
NPPmKIWCBEST 
POO.  5 

T5LO-  (NUM-1)  /2.0 
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T5HI«(N0M-l)/2.0 
CALL  WRITEO 
HUNT  ' 

REID  BN 
END 
C 

SUBROUTINE  POLYREOON ( ZERCR,  ZEROI,  G,  DEGREE,  K,  HJESTCNE, 
1  ®  ESnWO,  MDNEBEST,  MWCBEST) 

C 

C  SUBROUTINE  ID  FORM  Hi  AND  H2 
C 

DOUBLE  PRECISION  ZERCRU024)  ,  ZEROI  (1024) 

DOUBLE  PRECISION  ttNE(1024) ,  A3WO(1024) ,  B(3) 

DOUBLE  PRECIS  ION  CONE  (1024) ,  CIWO(1024) 

DOUBLE  PRECISION  HONE  (1024) ,  H3WO(1024)  ,  HMAX,  HMIN 
DOUBLE  PRECISION  ffiESTDNE (1024)  ,  ffiEST3WO(1024) 

DOUBLE  PRECISION  POM,  FOMQLD,  FOMDNE,  EOMWO,  AA,  BB 
DOUBLE  PRECISION  PRAT,  HBO0H(1024) 

INTEGER  G(1024) 

INTEGER  DECREE,  K,  ZCDUNT,  GCDUNT,  II,  JJ,  JJJ 
INTEGER  MDNE,  MWO,  MDNEBEST,  MWCBEST 
C 

C  OBTAIN  bl  AND  h2 
C 

MDN&O 

MEWOO 

DO  11  I>1,1025 
ACNE(D*O.ODO 
ATWO(L)"O.ODO 
11  CONTINUE 
C 

C  SELECT  A  ZERO 
C 

ZCDUMM. 

GOOUNIV1 

DO  5000  JJ-1,  EBGREE 
AApZERQR(JJ) 

B^ZEROI(JJ) 

C 

C  IF  DOG  PART  IS  VERY  SMALL,  LET  THIS  BE  A  LINEAR  EACTOR 
C 

IF  (DABS(BB)  .LT.1.0D-10)  BWJ.ODO 
C 

C  IF  IT  IS  THE  COMPLEX  CONJUGATE  ROOT  (I.E.  IMSG  PART 
NEGATIVE) 

C  THEN  SKIP  IT,  SINCE  IT  WILL  GET  PICKED  UP  BY  THE 
POSITIVE 

C  IMMG  QUADRATIC  FACTOR. 

C 

IF  (BB.LT.O.ODO)  GO  TO  5000 
C 

C  CREATE  A  QUADRATIC  FACTOR  FROM  A  COMPLEX  SET  OF  ROOTS, 


w  m  my 


C  A  LINEAR  FACTOR  FROM  A  REAL  ROOT. 

C 

B(l)-AA*AAfBB«BB 
IF  (BB.HQ.0.0D0)  B(l)— AA 
B(2)— 2.0D0*AA 
IF  (BB.BQ.0.0D0)  B(2)«1.0D0 
B(3)-1.0D0 

IF  (BB.  EQ.0.0D0)  B(3)«O.ODO 
C 

C  FOR  THE  SPECIAL  CASE  OF  TEE  ENTIRE  TRANSFER  FUNCTION  CN 
C  ONE  TRANSDUCER  ONLY,  LET  hi  BE  AN  IMPJLSE  (4) . 

C 

IF  (K.EQ.O)  TEEN 
MONEVl 

BCNE(l)  “1.0D0 

BID  IF 
C 

IF  (GCDUNr.GT.K)  GO  TO  232 
C 

C  INCORPORATE  THE  LINEAR  OR  QUADRATIC  FACTOR  INTO  THE 
C  hi  POLYNOMIAL 
C 

IF  (ZaXJNr.E3-G(GG0UNT))  THEN 
CALL 

POLMJLT(MONE,XWE,B,GONE,  ZCDDMT,  GQOUNT,  BB,  BONE) 
GCDUN1M3CDUN1M. 

GO  TO  5000 

BID  IF 
C 

C  INCORPORATE  THE  LINEAR  OR  QUAERATIC  FACTOR  INTO  THE 
C  h2  POLYNOMIAL 
C 

232  IF  (ZCODNT.NE.G(GCOUNT) )  THEN 
CALL 

POLYMOLTGONOr  ATWO,  B,  CIWO,  ZODUNT,  GODUNT,  BB,  HTWO) 

BID  IF 
C 

5000  Q0NT3NDE 

C 

C  SCALE  THE  COEFFIOBITS  TO  MAXIKJM  OF  1 
C 

HMAX-0.0D0 
DO  96  I^1,MDNE 

IF  (DABS(BCNE(U )  .GT.HMAX)  HMAX4ABS(BCNE(L) ) 

96  CONTINUE 

CO  97  LpI/MCNE 
HCNE(L)  HBCNE(L)  /BAX 

97  CONTINUE 
HMAXM3.0DO 

DO  98  1^1, mwo 

IF  (EABS(BIWOaJ)  .GT.HMAX)  HMAX*OABS(H3WO(U  ) 

98  CONTINUE 

00  99  I^l,fCWO 


wow 


fflWO(L)  =fflWO(L)  /HMAX 
99  CONTINUE 
C 

C  DETERMINE  THE  FIGURE  OF  MERIT  FOR  THE  DESIGN 
C 

IF  (K.H2.0)  THEN 

HCNE(l)  -l.ODO 
FOMXD°O.OBO 

END  IF 

DO  348  Iri.,MDNE 
ffiO!lH(L)-flDNE(U 

348  CONTINUE 

DO  349  M.+fCNE,  JENEH-MIWO 
mOIB  (L)  -HIWO(L-MDNE) 

349  CONTINUE 

CALL  HSTKE5(ffiOTB,M3NBfMIWO,FOM) 

C 

C  IF  IBIS  IS  THE  BEST  DESIGN  RELATIVE  TO  ALL  EAST  ONES, 

C  SAVE  OHE  RESULTS. 

C 

IF  (FOM.GT.FONQLD)  THEN 

HUNT  V**  BEST  YET  FOLLOWS  **' 

FQM3LDHFOM 

MONEBESTMENE 

KEWCBESTMOWO 

HUNT  *,,FO^,,FOM,‘  R»',K 

DO  111  II-1,ICNE 

ffiESTONE(II)  -BCNE(II) 

111  CONTINUE 

DO  2222  II«l,MWO 
ffiESnWO(II)  -BIHO(II) 

2222  CONTINUE 

C 

C  OPTIONAL  (DDE  TO  ERCVECE  SPECIFIC  INFORMATION  CONCERNING 
C  THE  NATURE  OF  EACH  SELECTED  ROOT  (I.E. ,  BASEBAND  OR 
C  SIDEBAND,  REAL  OR  CDMELEX) 

C 

c  ZOOUNT-l 

C  GCDUN3V1 

c  DO  1234  II-l, DEGREE 

c  AAfZERCR(II) 

c  HE^ZEROI(II) 

c  IF  (QABS(BB)  .LT.1.0D-10)  BB-O.ODO 

C  IF  (BB. LT.O.ODO)  GO  TO  1234 

c  IF  (GOOUNT.GT.K)  GO  TO  217 

c  IF  (ZOOUNT. H2.G(QOOUNT) )  THEN 

C  HUNT  Vffl.  ZERO: '  ,AA,  '  +/-',BB 

C  AA^DABS  ( 1 .  ODO-DSQPT  (AA*AAfBB*BB) ) 

C  IF  (AA.LT.1.0D-4)  HUNT  *,  'STOIBAND' 

c  IF  (AA.  GE.1  .OD-4)  HUNT  VPASSAND* 

C  IF  (BB.B3.0.0D0)  HUNT  V REAL' 

C  HUNT  V  ' 

c  ZCDUNP»ZOXJNPfl 
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c 
c 

C  END  IF 

c  217 
c 
c 
c 
c 
c 
c 
c 

C  END  IF 

c  1234 

END  IF 
C 

RETURN 
C 

END 
C 

SUBROUTINE  HSTATS(X,  M,  POM) 

C 

C  THIS  SUBROUTINE  DETERMINES  THE  STATISTICS  OF  THE  SAMPLES 
C  AND  RETURNS  THE  FIGURE  OF  MERIT  (FOM) 

C 

DOUBLE  PRECISION  X(1024) ,  XMAX,  XMIN,  XAVG,  XVAR 

DOUBLE  PRECIS  ION  POM/  XSUM/  XDUM,  XRAN3E 

INTEGER  M 

XSUFHJ.ODO 

XMAXm.ODO 

XMDKL.0D20 

XDUJHJ.ODO 

DO  10  1*1 /M 

XDU»=XDUMfX(I)*X(I) 

XSU»=XSUMfDABS(X(I)) 

IF  0QABS(X(D)  .GT.  XMAX)  XMAX*DABS(X(I) ) 

IF  OC(I)  .H2.0.0D0)  GO  TO  10 
IF  (DABSCX(I) )  .LT.XMIN)  XMUMJABSCXd) ) 

10  CONTINUE 

XAVOXSUM/EBLE  (M) 

XVAft*  (CBLE  (M)  *XDUM-XSUM*XSUM)  /  (EBLE  (M)  *  (EBLECM)  - 
l.ODO)) 

XRANjEVXMAX-XMIN 

IF  CXVAR.  EQ.O.ODO)  XVAi^l  .0D-10 

FOfrXAVQ'  (XRANGE*XVAR) 

35  RETURN 
END 
C 

C  POLYNOMIAL  RECONSTITUTION  SUBROUTINE.  THIS  ROUTINE  TAKES 
C  THE  LINEAR  OR  QUADRATIC  FACTOR  AND  MULTIPLIES  IT  BY  THE 
CURRENT 
C  POLYNOMIAL. 

C 


GCOUNPKjOOUNTfl 
GO  TO  1234 

IF  (ZOOUNT.  NE.G(GQOUNr) )  THEN 
PRINT  *,  *H2  ZERO:',AA/' 

AAfDABS  ( 1 .  ODO-DSQRT  (AA*AArfBB*BB) ) 

IF  (AA.LT.1  .OD-4)  PRINT  V  STOPBAND' 
IF  (AA.GE.1.0D-4)  PRINT  V  PASS  AND' 
IF  (BB.EQ.0.0D0)  PRINT  *, 'REAL' 
PRINT  *, '  ' 

ZCDUNT=ZaOCJNIH-l 

CONTINUE 


SUBROUTINE  FOLYttJLT  CM,  A,  B,  C,  ZCDUNT,  GdXJNT,  BB,  H) 


DOCBIE  IRBOSICN  A(1024)  ,  B(3)  ,H(1024) 

DOUBLE  IRBOSICN  C(1024) ,  BB 
INTEGER  Z COUNT,  H 
C 

C  IP  NO  CURRENT  POLYNOMIAL  YET  EXISTS,  1HEN  1HE  FACTOR 
BBOOMES 

C  THE  CURRENT  POLYNOMIAL. 

C 

IF  (M.B3.0)  THEN 

DO  1002  >1,3 
A(ti-B(L> 

C  (L)-B(L) 

1002  CONTINUE 

IML 

GO  TO  1421 

END  IF 
C 

C  IF  THE  CURRENT  POLYNOMIAL  IS  OF  ORDER  GREATER  THAN 
THREE,  THE 

C  RECURSIVE  RELATIONSHIP  APPLIES. 

C 

IF  (M.GT.3)  THEN 

DO  1001  >3,M 
C(L)-B(3)*A(L-2) 

C(U  -C(L)  4B(2)  *A(L-1) 

C(U  *C(L)  4B(1)  *A(U 
1001  CONTINUE 

END  IF 
C 

C  DUE  TO  VARIABIE  ADDRESSING  LIMITATIONS,  TAKE  CARE  OF  THE 
TWO 

C  HIGH  CREEK  AND  THREE  LCWEST  ORDER  COEFFICIENTS  MANUALLY 
C 

C(Mf2)-B(3)*A(M) 

C(MH)-B(3)  *A(M-1)  4B(2)  *A(M) 

C  (3>-B(3>  *A(1)+B(2)*A(2)4B(1)*A(3) 

C(2)  -B(2)  *A(1)  +B(1)  *A(2) 

C(D-B(1)*A(1) 

C 

C  UPDATE  THE  NUMBER  OF  ZEROS  USED. 

C 

1421  ZOOUNlVZOOUNIH-1 

C 

C  ORDER  INCREASES  BY  TWO  FOR  A  QUADRATIC  FACTOR,  CNE  EQR  A 
C  LINEAR  FACTOR. 

C 

IF  (BB.NE.0.0D0)  *"Mt2 
IF  (BB.BQ.0.0D0)  JMW-1 
C 

C  ESTABLISH  THE  NB9  CURRENT  POLYNOMIAL 


~f 
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A(D-C(U 
h(l)k:(U 
131  0GMT1NUE 
RETUro 
END 


w 
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